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ABSTRACT 


Four possible estimators are investigated for the monitoring 
of crustal deformations from a combination of repeated baseline length 
measurements and adopted geophysical models, particularly an absolute 
motion plate model. The first estimator is an extension of the familiar 
free adjustment. The next two are Bayesian type estimators, one weak 
and one strong. Finally, a weighted constraint estimator is presented. 

The properties of these four estimators are outlined and their physical 
interpretations discussed. 

A series of simulations are performed to test the four estimators 
and to determine whether or not to incorporate a plate model for the 
monitoring of deformations. It is concluded that it is preferable to 
adopt even a weak but realistic model than none at all. In this case, 
the weak Bayesian estimator (Best Linear Estimator — BLE) is preferred. 

It filters the signal (deformations) from the measurement noise in an 
optimal manner and, furthermore, is not overly sensitive to the errors 
in the adopted geophysical deformation model. 

The application of these estimations to the maintenance of a 
new conventional terrestrial reference system is discussed. The 
functions of the system are twofold. The first is to monitor the global 
rotations and translations of the earth with respect to an inertial frame. 
The second is to monitor the nonglobal motions or deformations of the 
earth. The relationship between these two functions is outlined. 
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1 INTRODUCTION 


"Length, Breadth and Thickness take up the whole of Space. 
Nor can Fansie Imagine how there should be a Fourth Local 
Dimension beyond these Three." 

(John Wallis, 1685) 


1.1 Background 

The relatively new concept of a dynamic and deformable earth 
and the recent availability of highly accurate geodetic measurement sys- 
tems lead us to address the problem of how to establish and maintain a 
suitable terrestrial reference system on such an earth. But first let 
us review briefly the history and status of the present system (see 
(Gulnot, 1978) for a more detailed account) and how it came to be con- 
sidered no longer satisfactory. 

In the latter part of the 19th century, the existence of polar 
motion was recognized. In order to monitor and correct for this 
phenomena which as it turned out could vary station coordinates by 
several meters, the International Latitude (ILS) was established in 
1899. The fundamental concept was to remove the systematic errors in 
the star positions from the latitude observations by choosing observa- 
tories on the same parallel (approximately 39“08') with identical 
instruments and procedures. Furthermore, it was assumed that the 
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stations were motionless with respect to each other and without varia- 
tions in their respective local plumblines. Thus, polar motion could 
be determined in a consistent and well defined manner by assigning con- 
ventional values for the initial latitudes of the five ILS stations. 

This convention was adopted in 1967 by the International Astronomical 
Union (lAU) and the International Union of Geodesy and Geophysics, 
defining the Conventional International Origin (CIO) for polar motion. 

The CIO pole definition, therefore, excluded the incorporation 
of observations from other stations. However, the polar motion deter- 
minations from only the ILS stations were not accurate enough particu- 
larly in monitoring short-term variations. This became especially 
apparent with the advent of new optical instruments and methods. 
Therefore, in 1962 the ILS was reorganized into the International Polar 
Motion Service (IPMS) which today provides pole coordinates at 0.05 year 
intervals based on latitude measurements from about 80 stations with a 
precision of about 1 meter. Here again a consistent set of initial 
latitudes is required and furthermore weighting procedures are applied 
since the quality of observations differ from station to station. 
Therefore, the IPMS pole differs from the CIO pole. 

The first axis of the present CTS is defined by the assigned 
longitudes of about 50 time observatories of the Bureau International 
de I'Heure (BIH) . The initial task of the BIH, created in 1912, was to 
maintain a uniform time scale. This function evolved, in addition, to 
monitoring variations in the earth’s rotation rate. Since this 
requires the pole position, the BIH began to compute its own polar 
motion values, hence a BIH pole. Furthermore, the BIH began applying 
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corrections to its computations based on earth rotation variations esti- 
mated from satellite Doppler and lunar laser ranging (LLR) observations 
(Feissel, 1980). These observations, as well as those from the other 
new space methods discussed below, are no longer referenced to the 
directions of the local plumblines as are the optical instrument obser- 
vations, but to terrestrial directions. 

Thus, today one has the choice of several sets of earth orien- 
tation parameters. The CIO pole adopted by the lAU and lUGG can no 
longer be considered accessible and of practical use (Kovalevsky and 
Mueller, 1981). The BIH provides the most frequent (five-day 
averages) and precise estimates (approximately OVOl (30 cm) for the pole 
components and 1 ms (45 cm) for UTl). However, today's requirements 
include estimates of polar motion every two days and variations in the 
earth's rotation rate each day, both with an accuracy of 5 cm or better 
(National Research Council, 1981). There is a further reason why the 
present situation is unsatisfactory. By coincidence, in the same year 
(1967) that the CIO pole was adopted, several scientific papers appeared 
in two distinct areas that effectively render one of the fundamental 
assumptions of the present system invalid. 

The theory of continental drift was first proposed by Wegener 
in 1912 although it was not accepted at the time by the scientific com- 
munity (Wegener, 1924; Uyeda, 1978). This idea was reborn fifty years 
later as plate tectonic theory. Culminating a previous decade of geo- 
physical investigations (McKenzie and Parker, 1967) and one year later 
(Morgan, 1968), advanced the hypothesis that the earth's lithosphere is 


3 


composed of rigid plates that are In relative motion along their 
boundaries, over an underlying as thenosphere. According to this theory, 
observations on the earth's surface are expected to have relative secu- 
lar motions on the order of 1-10 cm per year, depending on which plates 
the particular observatories are located. This, of course. Invalidates 
the CIO-BIH assumption of no relative motions between stations. 

This theory alone did not In Itself cause the lAU and lUGG to 
defer their CIO decision because at the time such small secular motions 
were lost within the noise of the available optical Instruments. How- 
ever, In that same eventful year results of the first experiments with 
Very Long Baseline Interferometry (VLBI) were published (Broten, et al, 
1967; Bare, et al, 1967). The opportunities presented by this measure- 
ment system cover nearly all aspects of reference frame requirements 
except the origin problem. First, recalling that radio Interferometry was 
Initially developed for astrometric applications, it has the potential for 
defining what is considered today, the most useful inertial reference sys- 
tem based on an adopted catalogue of extra-galactic radio source coordi- 
nates (Robertson, 1981). Second, it has the ability to measure interconti- 
nental baselines with near centimeter repeatability (e.g.. Herring, 
et al, 1981), and can be expected to monitor the relative motions 
between the observatories. Finally, it provides accurate estimates of 
short term variations in polar motion and earth rotation (e.g., 
Robertson, et al, 1978; Fanselow, et al, 1978). VLBI together with the 
other modern systems developed over the last 10-15 years. Satellite 
Laser Ranging (SLR) - (e.g.. Smith, et al, 1978), LLR (e.g., Mulholland, 
1978) and the Global Positioning System (GPS) - (e.g.. Fell, 1981) all 
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having the potential to meet the present earth orientation accuracy 
requirements, mark the urgency of defining a new reference system. 

The geodetic community reacted quickly to these developments in 
theory and instrumentation and since that time the problem of defining 
reference systems for a deformable earth has been the impetus for much 
activity. Two lAU colloqula (Ko^aczek and Weiffenbach, 1975; Gaposchkin 
and KoXaczek, 1981) have been entirely devoted to this topic and sum- 
marize developments in this area. A general consensus has emerged on 
how to remedy the present situation which will serve as a guideline 
for this investigation. 

It is generally agreed that two reference systems are needed. 

The first, using the nomenclature of (Mueller, 1981), is a Conventional 
Terrestrial System (CTS) in which positions and deformations on the earth 
could be described. The second is a Conventional Inertial System (CIS) 
to which the rotations and translations of the CTS could be referred. 

The first one is the focus of this investigation. 

The problem on the deformable earth is to establish an adequate 
representation of the earth's crust both spatially and temporally. We 
will follow the kinematic approach as outlined, for example, in 
(Gaposchkin, 1981; Kovalevsky and Mueller, 1981) which is the most con- 
ceptually simple and unambiguous. As in the CIO-BIH system, the CTS 
should be attached to observatories on the earth's surface. However, in 
the new system the stations cannot be assumed to be motionless with 
respect to each other. Furthermore, they should not be tied to the 
direction of the local plumblines but rather to directions tied to 
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the crust. These stations define a polyhedron whose edges or baseline 
lengths are accurately and directly estimable from VLSI and laser 
ranging observations. An adopted set of coordinates for these stations 
at an arbitrary initial epoch define what will be termed the fundamen- 
tal polyhedron. Implicit in the coordinates of its vertices are con- 
ventional spatial Cartesian axes that define the reference frame. These 
are accessible at any epoch through global transformations (rotations 
between the true frame to which the nutations refer^ and translations 
with respect to the initial origin) . The problem is to relate the 
rotated, translated and deformed polyhedron at a later epoch to the 
fundamental polyhedron, or equivalently to both the CIS and CTS. 
Therefore, the functions of the CTS stations are twofold. First, an 
extension of the present BIH system is to monitor those motions common 
to all points on earth (the polyhedron - CIS connection) . The second 
function is to monitor what is left, i.e., the deformations of the poly- 
hedron (the polyhedron - CTS connection). Therefore, by definition, the 
deformations do not contain any common rotations or translations (that 
are statistically significant) . 

It is conceptually possible to directly monitor the orientation 
parameters between the CTS and CIS frames simply by 3 rotation angles 
between two purely kinematlcal reference frames (the CIS frame being 
defined by a unit sphere polyhedron of extra-galactic radio sources - 
Gaposchkin, 1981) . These rotations include the combined effects of 
precession, nutation, polar motion and earth rotation. However, since 
precession and nutation can be computed quite accurately from adopted 
earth models, it is preferable to model these effects. The earth 
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orientation parameters (polar motion and earth rotation) which at this 
stage cannot be modeled adequately are estimated through observations. 

Of course, any errors in precession and nutation as well as common rota- 
tions due to plate motions will be absorbed in these parameters. The 
latter could occur from an inadequate distribution of observing stations 
(for example, all stations on one plate). Considering that the require- 
ments for earth orientation parameters call for 1-2 day resolution at 
the subdecimeter level, near continuous observations are required. How- 
ever, since these parameters are global, it is possible that only a 
subset of the CTS stations will have to participate on a regular basis. 

The orientation of the CTS reference frame axes are quite arbi- 
trary. However, it is widely agreed that efforts should be made to 
maintain continuity with the BIH system. This could be accomplished 
most simply through VLBI observations as will be described in Chapter 2. 
Similarly, the origin could be defined arbitrarily but it is preferable 
that it be at the center of mass of the earth, as determined by satel- 
lite or lunar laser observations. Changes in the center of mass could 
be determined in the same way, or by absolute gravity measurements 
(Heiskanen and Moritz, 1967; Mather, 1973; Moritz, 1979; Zielinski, 
1981). 

For the monitoring of deformations, it is anticipated that dis- 
placements due to the tidal potential and loading effects can be com- 
puted to within centimeters and therefore will be modeled (Mueller, 
1981). Furthermore, although the stations should be located on stable 
regions of their respective plates, any possible site stability prob- 
lems should be monitored by on-site instruments such as gravimeters and 
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by local geodetic networks, particularly using GPS interferometry 
(Counselman and Shapiro, 1978). The deformation monitoring functions 
of the CTS stations, then, will primarily be directed to interplate 
motions. As for modeling plate tectonic motion, this is a controver- 
sial question since it is uncertain whether the average long-term plate 
motions Inferred from geophysical and geological data relate to current 
rates of motion. The adoption of a plate model has been advocated most 
strongly in (Bender, 1974; Bender and Goad, 1979; Bender, 1981) and 
will be addressed quite extensively in this investigation. 

The fundamental polyhedron as defined by the adopted CTS station 
coordinates at an initial epoch has a certain size and shape. At a 
later epoch, the deformed configuration is completely determined from 
the changes in the baseline lengths which are thus the key to monitoring 
def onnations . Since the CTS at any epoch is defined by the coordinates 
of all the stations, and considering that deformations due to plate 
motions are predicted at the decimeter per year level, only periodic 
re-observations of the baseline lengths need to be taken but from all 
stations (including the ones that monitor earth orientation on a regu- 
lar basis) . 

The realization of the CTS should have low sensitivity to 
changes in the distribution of the observing stations and in the fre- 
quency of observations from individual stations, considering that the 
number of stations and observations is likely to change from time to 
time. It should avoid as much as possible any dependence on geophysi- 
cal hypotheses although it is recognized that the CTS will require some 
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geophysical information, for example (earth tide or plate motion para- 
meters as mentioned above). The definition of the CTS, however, should 
not depend on this information. 

The definition of the reference system should be compatible with 
simple operational descriptions of how the system should be utilized. 
There should be established procedures or algorithms for acquisition, 
reduction, and application of observational data. This includes the 
adoption of earth models and fundamental constants to be adhered to in 
all computations. 

1.2 Purpose and Ob.lectives 

It now appears that a new CTS (and CIS) will soon be estab- 
lished possibly (and hopefully) as an outcome of the upcoming 1983 
MERIT (to Monitor E^arth Rotation and ^ntercompare the techniques of 
observation and analysis) main campaign (Wilkins, 1981). In anticipa- 
tion of this event, this investigation addresses certain aspects of the 
problem of how to define and maintain a new CTS on the deformable earth. 

The approach to be followed in this investigation, using as 
guidelines the consensus described in the previous section, is outlined 
in Fig. 1. The fundamental polyhedron and the reference frame are 
defined at an initial epoch t^ by the adopted (fundamental) coordinates 
X of a particular set of stations. The coordinates are to be esti- 
mated from a dedicated observation campaign of VLBI, SLR, LLR and pos- 
sibly GPS stations as described in Chapter 2. 

The level I stations include a dedicated subset of the CTS 
observatories which monitor earth orientation parameters (EOP) on a 
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10 











continuous basis. The polar motion and earth rotation parameters are 
averages over say 1-2 day periods as indicated by the short intervals in 
Fig. 1. They are rotations between the true frame to which the nuta- 
tions refer and the reference frame axes defined by X . In order to 
insure that these parameters continue to refer to the same set of axes, 
the deformations of the polyhedron need to be estimated periodically but 
from measurements at all the CTS stations (level II) as indicated by 
the shaded portions of Fig. 1. For example, at the end of observation 
interval A, the deformations AX are estimated as will be described in 
Chapter 3 and are then added to the fundamental coordinates X to 
define the CTS for the next observational period B . In this way, the 
ref erence system is maintained, i.e., the earth orientation parameters 
refer to the same reference frame defined by X . As indicated by the 
addition of the terms in brackets to X , we allow the possibility of 
updating the station coordinates on the basis of an adopted plate motion 
model. This could also be viewed as a correction to be added retroactively 
to refine the earth orientation parameters estimated over Interval A. It 
should be mentioned that the number and distribution of the CTS stations 
of level I can vary without affecting the reference system, assuming 
that the deformations are estimated from a well distributed network of 
level II CTS stations. 

As mentioned in the previous section, any changes in the earth’s 
center of mass, relative to that defined by X , can be determined from 
satellite dynamics or absolute gravity measurements. In addition, a 
scale parameter which would represent an expansion or contraction of the 
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earth could be determined from re-measured baseline lengths or from 
gravity observations (Helskanen and Moritz, 1967). Both the transla- 
tional and scale parameters unlike the rotational ones, would have to be 
monitored from a well distributed global network, preferably from all 
CTS stations . 

It should be clear that the reference frame defined by X , once 
chosen, does not change. It can be thought of as fixed to the initial 
(at tg) positions of the CTS stations (the fundamental polyhedron) . The 
frame, therefore, consists only of a set of spatial Cartesian axes with a 
particular orientation and origin. It is the reference system that is 
changing and moving with the deformed polyhedron. The fundamental poly- 
hedron-CIS connection is given by the EOF. The fundamental 
polyhedron-CTS (deformed polyhedron) connection is given by the esti- 
mated deformations. Therefore, the CTS and its frame coincide in 
general only at the initial epoch. However, the CTS is not only a set 
of changing station coordinates. It must contain a well-defined 
description of anything that would influence these coordinates. This 
includes, of course, X (the reference frame), the CTS stations, adopted 
earth models (precession, nutation, tides, plate tectonics, etc.) and 
related fundamental constants, parameter estimation models and estab- 
lished procedures for all CTS operations as mentioned in section 1.1. 

The purpose of the CTS, then, is to make the reference frame accessible 
to the user who can then determine positions and detect motions of the 
earth. 
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It is useful to present here an excerpt from the concluding com- 
ments from the 56th lAU Colloquium (Gaposchkin and Kolaczek, 1981) as 
summarized by Kovalevsky and Mueller. These points concern the actions 
required before final decisions are made with regard to the new CTS. 

They include 

1. Selection of observatories whose catalogue will define the CTS. 

2. Initiation of measurements at the observatories. 

3. Recommendation on the observational and computational mainte- 
nance of the CTS (e.g. permanent versus temporary and repeated 
station occupations , constants to be used) . 

4. Decision on how far and which way the earth deformation should 
be modeled initially. 

5. Plans and recommendation for the establishment of new interna- 
tional services to provide users with the appropriate informa- 
tion regarding the use of the CTS frame. 

The objectives of this investigation address certain aspects of points 
1, 3 and 4 listed above. The primary focus will be on the choice of 
algorithms for the estimation of earth deformations and how this relates 
to the estimation of earth orientation parameters and the maintenance of 
the CTS (point 3) . The choice of algorithm will depend on the question 
of whether or not plate motion models will be adopted. An attempt is 
made to answer the question of whether adding a model will improve the 
deformation monitoring capabilities of the system (point 4). This will 
in turn Influence somewhat, as will be investigated, the selection of 
observatories whose coordinates will define the CTS (point 1) . 


1.3 Organization and Scope 

Chapter 2 discusses the definition and maintenance of the CTS. 

A method Is outlined for estimating the coordinates of the fundamental 
polyhedron from a network containing different measurement systems. 

Next, the separation of global and deformation parameters Is addressed, 
that Is the polyhedron-CIS and polyhedron-CTS connections. An example 
with a VLBI network Is given for which a particularly suitable para- 
meterization for CTS operations Is proposed. Different approaches are 
discussed for maintaining the reference frame. A set of constraints 
are derived that maintain a discrete Tlsserand's mean axes of crust. 

The physical Implications of adopting an absolute plate motion model 
are described. 

In Chapter 3, four estimation algorithms are presented for the 
analyses of polyhedron deformations. Two are chosen from the class of 
biased estimators considering the Inherent singularity In maintaining a 
reference frame solely from geodetic observations. Two conditionally 
unbiased estimates are then presented. The optimal properties of these 
four estimators are outlined as well as their physical meaning. All are 
able to deal with the Incorporation of geophysical data, particularly 
provided by an adopted plate tectonic model. In fact, three of the 
estimators require some a priori deformation Information In addition to 
the geodetic observations. Finally, least squares collocation Is 
applied to deal with the addition or temporary loss of several CTS sta- 
tions . 
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In Chapter 4, numerical tests are performed to gain a better 
appreciation of the optimal properties of each estimator. The main 
thrust of these tests is to determine the best deformation estimation 
algorithm, of the four presented in Chapter 3, for combining geodetic 
and geophysical data. In anticipation of the planned MERIT main cam- 
paign scheduled for 1983-84, we extend the optimal polyhedron design 
analysis of (Mueller, et al, 1982). There it was assumed that no geo- 
physical model for deformations is available. Here, we study the 
effects of adopting an absolute motion plate model. 

The appendices outline some recent results in the approximate 
theory of optimal design which are applied to the polyhedron design 
problem. In addition, several results are presented dealing with 
weighted pseudoinverses. 


15 


2 DEFINITION AND MAINTENANCE OF THE CTS 


2.1 Introduction 

In this chapter, a method is outlined to separate global and 
deformation parameters in the maintenance of the CTS. Recall that 
it is necessary to monitor the deformations and update the initial coor- 
dinates of the CTS stations so that the transformation parameters (glo- 
bal rotations and translations) will always refer to the same reference 
frame. In order to do this, either some constraints are needed to 
insure that the deformations indeed do not contain any global motions, 
or some external information on the expected deformations is required. 
Both of these approaches are discussed. In the former, the reference 
frame axes are seen to be a discrete version of a Tisserand*s mean axes 
of crust. In the latter approach, the reference frame axes are fixed in 
that part of the earth to which the a priori deformations refer, for 
example the mantle when using absolute plate motion models. We begin 
with a method to estimate the fundamental coordinates (we drop the 
subscript t of Chapter 1) that define the reference frame. 

2.2 Definition of the CTS Frame 

As discussed in Chapter 1, the main need for a new CTS comes from 
the anticipated accuracy of the new geodetic measurement systems and 


16 


subsequently their ability to monitor deformations of the magnitude pre- 
dicted by plate tectonic theory. The systems can be divided into base- 
line methods and coordinate methods. In the former, primarily 
interferometric methods such as VLBI and GPS interferometry, the origin 
of a terrestrial coordinate system is not sensed by the observables. 

The smallest unit for this type of measurement is a baseline (hence, 
baseline method) although it is a triangle for earth orientation analy- 
sis (Bock, 1980). The estimable parameters are baseline lengths and 
variations in earth orientation (relative to a well-defined initial 
orientation, i.e. a CTS frame). 

The coordinate methods include SLR, LLR and GPS in the Doppler 
and pseudo-ranging modes (e.g.. Fell, 1980). The basic unit for these 
systems is one station. An origin is sensed and thus these systems 
provide Cartesian coordinates. SLR observations are sensitive to polar 
motion variations but basically insensitive to variations in earth rota- 
tion (Van Gelder, 1978). On the other hand, LLR observations are sensi- 
tive to all three earth orientation parameters but particularly to 
variations in earth rotation (Calame, 1982). 

In all of the new systems, though, coordinates (or in the case 
of interferometry, coordinate differences) are inseparable from earth 
orientation parameters. This means that to estimate coordinates (or 
coordinate differences) an external source of earth orientation is 
required and vice versa. By using BIH earth orientation values the CTS 
frame can be made continuous with the present BIH frame at the initial 
epoch as will be described below. 
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Considering the different measurement systems available, it will 
be necessary to merge several networks, each one defining essentially its 
own reference frames both CTS and CIS, into a common set. Suppose the 
relations between two CIS's is 


= R^(aj^)R2(a2)R3(a2)x^ 

Similarly, the relation between two CTS's is 
= R3^(8^)S2(32)R3(e3)X^ 

The transformation from CIS to CTS is (Mueller, 1969) 



S^NPx^ 
S NPx 


II 


( 2 . 2 - 1 ) 


( 2 . 2 - 2 ) 


(2.2-3) 


where common nutation (N) and precession (P) matrices are assumed to be 
used in both techniques. The earth orientation matrix is given by 


S = R 2 (-C)Ri(-n)R 3 ( 6 ) (2.2-4) 

in which ^,ti are the coordinates of the pole and 0 is the Greenwich 
Apparent Sidereal Time. After some reduction and neglecting 
second-order terms, we have from (Mueller, et al, 1982) 

-An = -(ti^ - n^^) = -6-j^ + aj^cos 0 + q 2 sin 0 (2.2-5) 

-AC = -(C^’-C^’’) = -82 " ttj^sine + a2cos0 (2.2-6) 

WAUTl = W(UT1^ - UTl^^) = - 83 + (2.2-7) 
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where W is the ratio of universal (UTl) to sidereal time. By station 
collocation, i.e., maintaining different instrument types at common 
sites, one determines the CTS difference (3 angles). Then through the 
earth rotation parameter differences one finds the CIS differences (a 
angles). This indirect approach has been suggested by (Kovalevsky, 
1980). The determination of the CIS differences is treated in (Mueller, 
et al, 1982). 

In the following, a method to estimate the CTS differences and 
simultaneously to define the CTS frame based on the estimation of a 
unique well-defined set of coordinates for the fundamental polyhedron is 
outlined. Suppose that one baseline method, VLBI, and two coordinate 
methods, SLR and LLR, are to participate in a campaign to estimate the 
fundamental global spatial Cartesian coordinates, X^. VLBI observations 
are Insensitive to the initial orientation of the CTS frame, making 
coordinate differences and earth orientation parameters (polar motion 
and UTl) Inseparable. In practice, this dependency is broken by ini- 
tializing earth orientation over, say, the first day of observations of 
a particular campaign. As shown, e.g., in (Bock, 1980), the estimation 
of baseline components is then biased by this initialization but in this 
case the bias can be used to our advantage. Continuity with the pres- 
ent terrestrial system can be achieved by the input of BIH polar motion 
and earth rotation values for the initial step. In this way, at the 
fundamental epoch t^ the new CTS frame can be aligned with the BIH frame 
through the estimated VLBI coordinate differences (within the errors in 
the BIH values) . The first axis (x) of the BIH frame is in the direc- 
tion of the Greenwich mean astronomic meridian and the third axis (z) is 
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towards the average (mean) north terrestrial pole. The second axis (y) 
completes a right-handed global spatial Cartesian coordinate system. 
From the SLR or LLR estimated coordinates, the origin of the CTS frame 
could be made geocentric. 

One is then led to the following three sets of transformation 
equations from which Xq can be estimated 


0 63 -63 63 

-63 0 63 (X3>3+ (2.2-8) 

[S2 -h U- ' 


■ ° ^3 -V 

(l+<^>(Xi)j+ -Y 3 0 Yi 

LT^2 “ - 

“V - (2.2-10) 

The first set (2.2-8), one for each LLR station, has as observations 

i 

the geocentric coordinates of site i . The parameters are the CTS 
LLR site coordinates three rotation angles ^3 (connecting 

LLR to VLSI) , a scale factor c^^ (LLR to VLBI) and three translation 
parameters 6 ^, 62 , 6 ^ to SLR origin). For the second set (2.2-9), 

the parameters are the CTS SLR site coordinates (X^)g» three rotational 
angles y^, Y 2 , y^ (SLR to VLBI) and a scale factor C 2 (SLR to VLBI). 

For the third set (2.2-10), the parameters are the VLBI site coordi- 
nates (X^)^, and the observations are any Independent subset of coordi- 
nates differences AX^ from the VLBI estimated parameters. Different 
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combinations of equations (2.2-8) to (2.2-10) could be formulated 
although those given here reflect today's situation (origin defined by 
SLR, orientation and scale by VLBI) . Additional sets of equations could 
be added for other measurement systems (e.g., GPS). Considering (2.2-8) 
to (2.2-10) as observation equations and computing a weight matrix from 
the covariance matrices of the laser ranging and VLBI adjustments, one 
could then perform a least squares adjustment to estimate X^, a consist- 
ent set of CTS coordinates (at the collocated sites ~ ^^i^S ~ 

(X^)^ is constrained) that would define the new reference frame. It 
would be geocentric and aligned with the BIH frame at t^. 

2.3 Separation of Global and Deformation Parameters 

As mentioned in the previous section, all the available modem 
measurement systems suffer from the position-orientation inseparability 
problem which is inherent in all strictly earth based observations. 
However, it is required for the definition of the CTS to estimate peri- 
odically station coordinate changes (deformations) which are free from 
global rotations and translations. In order to separate global motions 
and deformations one approach is the addition of a set of constraints to 
overcome the singularity problem. These constraints are derived in the 
next section. 

Without any external information of expected deformations, one 
is led to a free adjustment with orientation parameters of which several 
approaches are possible (Fritsch and Schaffrin, 1981) . Here a two step 
analog to the "classical approach" is presented. First, the earth 
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orientation parameters are eliminated, more precisely those motions com- 
mon to all stations. Then, the deformations are estimated from periodi- 
cally repeated measurements of the baseline lengths by minimizing a 
weighted norm in the parameter space of coordinate changes. Simul- 
taneous adjustments (one step approaches) have been proposed for single 
measurement types particularly VLBI (Cannon, 1979; Manabe, 1982). An 
example, with a VLBI network will help elucidate these ideas. 

Consider the basic VLBI mathematical model as given in (Bock, 
1980) . For baseline i-j observing source 1 at epoch k the path dif- 
ference (time delay times the speed of light) which the Incoming signal 
must travel after its reception at station i till its arrival at 
station j is given by 


ijk£ 
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+ c(ACf> + Ac, (t -t^)) 

Ui K u 


(2.3-1) 


where (a, 6) are the source coordinates, ACq,Ac^ are clock parameters 
and c the speed of light. As mentioned above, the coordinate dif- 
ferences Ax, Ay, Az and the earth orientation parameters of polar 
motion (C,ri) and earth rotation (6) are inseparable. This can be seen 
by an examination of the linear relationships between the coefficients 
of the design matrix A 
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(2.3-2) 


A. = 


AxA, - AzA. 
Az Ax 


A = -AyA. + AzA. 

Ti Az Ay 


A- = -AyA. + AxA, 

0 ■'Ax Ay 


where the A’s denote the partial derivative of the observable with 

respect to the subscripted parameter which can be found in (Bock, 1980). 

There, a different parameterization is proposed in which the total 

observation period is split up into earth orientation steps. For the 

m'th step, 3 rotation angles ^ , H - ri, » 6 “6, are estimated, 

m 1 m 1 m 1 

relative to one fixed step, where 6^ are input into the adjust- 

ment as determined from external sources (e.g. ,BIH values). That is, 
the earth orientation parameters are averages over certain time inter- 
vals (e.g., 1 day averages). One set of coordinate differences are also 
estimated for the entire observational period. This parameterization is 
useful in some situations, one was mentioned in the previous section. 
However, it is particularly unsulted for CTS operations since the func- 
tion of the CTS is to separate earth rotation variations and station 
displacements in a well-defined manner. In this parameterization, it is 
difficult to maintain orientation continuity and furthermore the size 
and distribution of observations of the first step are arbitrary. More- 
over, this fixed step biases the estimation of coordinate differences by 
any errors dri^^, dG^^ in the external earth orientation information 

according to 
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(2.3-3) 
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where T, e, a are the "biased" baseline components (Bock, 1980) 

The one step free adjustment approach for VLBI is as follows. 
Consider the parameter vector 



(2.3-4) 


where represent the coordinate changes , X 2 the earth orientation 
parameters (3 per step) and X^ other "non-geodetic" parameters such as 
the right ascension differences (a^-a^) , declinations (6^) and clock 
parameters (Cq,Cj^) appearing in (2.3-1). The normal equations can be 
written as 
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(2.3-5) 


where 



A^PAj 


(2.3-6) 




(2.3-7) 


P is the weight matrix of the observations L . These are rank defi- 
cient by six due to the origin (VLBI is insensitive to an origin) and 
orientation (coordinate-orientation inseparability problem) defects. 

One proceeds in '^two steps. First, the orientation and model 
parameters are eliminated such that 
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■==2] p*22 r”2‘ 

-*3. -*'23 ’*33- -®3- 


(2.3-8) 


Note that only the Cayley inverse is required in this step. Next, and 
most generally, a weighted norm is minimized 

II .. II /-..^TLfVN ^ /O O 


xllj^ = (X MX)' 


(2.3-9) 


where M is a weight matrix constructed from an adopted deformation 
model as described in Chapters 3 and 4. Denoting 


N = Nil - [N^2 


^22 «23 


”23 *^33 


U = Ui - [N^2 N 13 I 


^22 "23r'K 


"L "3 


the coordinates of the CTS stations are estimated by (see section 
3. 3. 1.2) 


= m"^n(nm“^n)'’’u 


where + denotes the pseudoinverse. Here the matrix M is assumed to 
be positive definite, the more general case of a positive semidefinite 
M is treated in Chapter 3. If there is no a priori information model, 
i.e. , M = I 


/V 

= N U 
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the ordinary pseudoinverse solution. This choice of M-norm for is 
equivalent (see the next section) to defining a discrete Tlsserand's 
mean axes of crust. Instead of using the pseudoinverse to compute 
a set of constraints derived and explained in the next section could be 
augmented to the set of normal equations. In the second step, X 2 and 
X^ are estimated by (2.3-8). 

As mentioned above, several other norm choices are possible. 
Some of these including the "classical" approach presented here are 
applied in (Dermanis, 1977, 1981; Manabe, 1982) although In the less 
general case of M = I. Other approaches are possible, too. One might 
consider X^ as a random variable vector and X 2 and X^ deterministic. 
This would lead to a least squares collocation approach (section 3.6). 
Other combinations are possible. In order to circumvent choosing from 
all these different possibilities, the number of which indicates the 
problems encountered from the coordinate-orientation inseparability 
problem, an analogous but unambiguous approach is outlined below. It 
takes into consideration that the X^ parameters need to be monitored 
periodically from all CTS stations while the X 2 parameters must be 
monitored continuously from possibly a subset of the stations. Fur- 
thermore, the impractical simultaneous adjustment of several observa- 
tion types (e.g., VLBI, SLR, LLR) is avoided. 

Consider the following two step approach continuing the VLBI 
network example. Once the fundamental coordinates X^ are adopted, the 
following parameterization is particularly suitable for the first step 
of monitoring earth orientation. Instead of parameterizing the coordi- 
nate differences that appear in the mathematical model of equation 
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(2.3-1), we parameterize baseline lengths for the set X^. The set X 2 
includes 3 orientation parameters per step (e.g., one day earth orien- 
tation and two day polar motion averages) , but in this case there is no 
need to fix one step. These earth orientation estimates are with 
respect to the defined reference frame axes implicit in the CTS coordi- 
nates that are input into the adjustment. Thus, continuity in orienta- 
tion is maintained consistently even with a gap in the observations and 
no bias is introduced as in the previous parameterization. 

The partial derivatives for the baseline lengths (the X^ set) 
can be derived from the coordinate difference partlals 


Ax 


ij 


Ay 


ij 


Az 


ij 


= -COsd^^COsOj^- 
= cos6^sin(6j^ - a^^) 
- -slnSj, 


Then, 


“ij 


Axfj 


(2.3-14) 


(2.3-15) 


The baseline vector can be written as 
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where a and 3 are baseline direction angles, from which 


(2.3-16) 
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(2.3-17) 


so that 
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(2.3-19) 


“ij 


[A Ax., + A. Ay + A. Az ]/£. 
Ax^^ ±2 •'ij Az^^ i ij 


(2.3-20) 


All the other elements of the design matrix A are the same as in the 
coordinate difference adjustment (Bock, 1980) . 

This parameterization allows a clean separation of rotation and 
deformations. The baseline lengths (X^) can detect any possible 
short term deformations but recall that only a subset of the CTS sta- 
tions is observing. Although some have argued that VLBI observations 
are sufficient for monitoring earth orientation (e.g., Gaposchkin, 

1981) other measurement types will most likely be involved, too. In 
that case, the CTS coordinates input into, say, SLR observation adjust- 
ments, define the reference frame to which the earth orientation para- 
meters refer. Recall that a consistent set of polyhedron coordinates 
at the initial epoch would be determined from several measurement sys- 
tems as outlined in section 2.2. These coordinates are updated 
periodically using the deformations estimated from the second step of 
this approach (see below) , in order that the earth orientation 
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parameters of the first step refer always to the same reference frame. 
Thus, the coordinates of the stations are never included as parameters 
for the earth orientation adjustments of any participating measurement 
system. In this baseline length-orientation angle parameterization, 
the orientation angles include any common rotations of the crust due to 
plate motions relative to the mantle (although this is not expected as 
discussed below) . The baseline lengths may include a global expansion 
or contraction. These types of motion will be absorbed in the CTS by 
definition (Kovalevsky and Mueller, 1981). 

For the periodic monitoring of deformations, the parameters of 
interest are the baseline lengths (X^) of the deformed polyhedron 
defined by all the CTS stations with possibly different measurement 
systems. The baseline lengths contain only deformation information 
since they are independent of coordinate system (assuming that the 
coordinate system is not changing over the interval of time in which 
the baselines are estimated, for example, over a long arc laser ranging 
solution) . They indicate the change in size and shape of the polyhedron 
relative to the fundamental polyhedron. Furthermore, the baseline 
lengths are the only strictly estimable parameters from all the dif- 
ferent measurement systems and also the most accurate ones. The proce- 
dure, then, is to pool the baseline length estimates from the separate 
adjustments of the respective measurement systems, along with their 
covariance matrices. In a second adjustment, which is the subject of 
Chapter 3, the deformations of the CTS stations can be estimated by com- 
paring the deformed baseline lengths to their initial epoch values. 
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It Is In the second adjustment that a plate motion model could 
be Introduced to improve the estimation of the deformations, and to test 
its consistency with the geodetic observations. If no plate motion 
model is adopted, then In order to maintain the CTS, a set of constraints 
is required in the estimation of deformations to Insure that they 
Include no global motions. These are derived in the next section and 
are generalized to Include geophysical models of deformations. In 
the presence of such models, other approaches to maintain the CTS will 
also be described. 

2.4 Approaches to Maintaining the CTS 

2.4.1 Reference Frame Constraints - Geometrical Approach 

The reference frame axes at an initial epoch are conventionally 
defined by the adopted coordinates of the CTS stations (the fundamental 
polyhedron). At a later epoch, the deformed polyhedron will similarly 
define a different set of axes. In order to maintain a consistent 
reference frame one approach is that these two sets of axes are con- 
strained to differ by only global transformations as least in some 
weighted least squares sense. In this section, the constraints for 
this approach are derived from geometric considerations. 

At the fundamental epoch the polyhedron coordinates are given 
by the initial coordinates X^, At a later epoch, the new coordinates 
from which global rotations and translations have been either modeled 
(precession and nutation) or estimated (polar motion and earth rotation) 
and removed, differ from Xq. Furthermore, any translations due to 
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changes in the center of mass have been corected for. Any remaining 
differences between these two sets of coordinates are due to deforma- 
tions. The constraints necessary to insure that the reference frame 
axes implicit in and are equivalent can be derived as follows 
(expanding on an example in (Leick and Taylor, 1980)). Consider that 
the two sets of coordinates for station i 



The scale problem will be addressed at a later stage. Here the assump- 
tion is that any scale differences between the various systems have been 
determined and reconciled as in section 2.2. The above equations can be 
rewritten as 3 observations equations per station 



where V is the noise vector. In matrix form. 
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L = AX + V 


(2.4-4) 


The "observation vector", L Includes the changes in the coordinates 
with respect to the initial epoch, i.e., the deformations. The "para- 
meter vector", X Includes the common rotation and translation com- 
ponents (what is left after removing precession, nutation, center of 
mass changes, etc.). Introducing a weight matrix P for the observa- 
tions we arrive at the standard least squares estimate (Uotila, 1967) 

X = (A^PA)~^a'’'^PL (2.4-5) 

In order for the reference frame implied in X^ and X^^ to be maintained, 
it is required that X = 0. This occurs in two non-trivial cases. 

First, when L = 0, which would occur only for a rigid earth (v/lthin the 
observational noise). The second, more interesting case would be when 

A*^PL = 0 (2.4-6) 

Now, let us reconsider the problem. In reality, the deforma- 
tions are the parameters to be estimated so that (2.4-6) becomes our 
reference system maintenance constraint 

CMX = 0 (2.4-7) 

Consider M as a weight matrix derived from an adopted model for sta- 
tion deformations as described in section 4.2. It defines a weighted 
norm in the parameter space 

II X lljj = (X'^MX)^/^ (2.4-8) 
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The former design matrix A is now the constraint matrix C which has 
the following form 


where 


C = 


S, S„ ... S 

12 p 

I I ... I 



(2.4-9) 


(2-4-10) 


I is the 3x3 identity matrix and p is the number of polyhedron sta- 
tions. To understand the significance of (2.4-7), let us assume for a 
moment that M is diagonal and the weight for the i'th station is m^. 
The constraints can then be divided into two sets 


0 


P 

S m, 
i=l 





(2.4-11) 


(2.4-12) 


These constrain, respectively, the orientation and origin defined by Xq 
and Xj to coincide, in the weighted least squares sense. 

Consider again the original problem. We know that the standard 
least squares estimate has the property that 


V^PV = minimum 
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(2.4-13) 


/V 

Since X = 0, it follows from 


A.»p /\ »T 

V PV = L PL - 
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(2.4-14) 


that 


L PL = minimum 


In our reversed problem this implies that 


ys, 

X MX = minimum 


(2.4-15) 


(2.4-16) 


i.e., a minimum norm solution. 

Thus, in the general case, we can say at this point that the 
reference frame axes are fixed in the earth's crust, through a discrete 
number of CTS stations, in a minimum M-norm sense. In the next chapter, 
we present estimation techniques that incorporate these constraints. 

The constraints (2.4-7) are given above in terms of global 
spatial coordinates. Let us express them in a local geodetic system 
(u,v,w) where the u-axis is positive north along the geodetic meridian, 
the v-axis is positive east and w-axis in the direction of the ellipsoi- 
dal normal. From (Rapp, 1976), the relation between the two systems is 
given by 
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Az 


(2.4-17) 


where <})^, are the geodetic coordinates of station i . Assuming 
again that M is a diagonal matrix and approximating the ellipsoid by 
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a sphere of mean earth radius, we arrive at the two sets of constraints 
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(2.A-18) 


(2.4-19) 


that correspond to (2.4—11) and (2.4-12), respectively. The constraints 
(2.4-18) correspond to those given in (Bender and Goad, 1979; Bender, 
1981) . The w coordinate (which .indicates geometric height) is in 
parentheses since it appears in the matrix multiplication. This means 
thht an infinitesimal rotation only causes horizontal coordinate changes , 
i.e. , the height is insensitive to such rotations. However, the 
opposite is not correct. Rotation is sensitive to height changes, the 
proof of which is due to an unpublished manuscript of S. Y. Zhu and is 
reproduced here. This will indicate that the constraints (2.4-18) and 
(2.4-19) should be taken as a complete set. 

An infinitesimal rotation a and a translation 6 will change the 
coordinates according to (2.4-3). The corresponding height change for 
a particular station is 


Aw = 


Ax * X + Ay « y + Az 

. 2 ^ 2 ^ 2 - 1/2 
[x +y +z J 


(2.4-20) 


If only an infinitesimal rotation exists, that is 6 = 0, then we will 
find Aw = 0 as stated above. 
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On the other hand. It can be shown that a height change Aw can 
influence rotation significantly. Consider (2.4-3) again as observa- 
tion equations where the parameters are the infinitesimal rotations a 
and translations 6, and the observations are changes in the coordinates. 
The normal equations are 


NX = U 


(2.4-21) 
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Compare the expression for U to (2.4-7) and (2.4-9). In general, as 
seen in (2.4-22), the off-diagonal elements of N are not equal to 
zero. We write N ^ as 
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Transforming the above equation into the local geodetic system and 
setting u = v=0 (see (2.4-19)), that is considering only a height 
change , 
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£Aw cos(|>cosX 
ZAw cos(t>slnX 
ZAw sin4> 


(2.4-24) 


and 


“l = *16“6 

Since a^^, a^^, a^^ of N ^ generally are not equal to zero, neither is 
a^. Simulation shows that can be of the same order ot magnitude as 
Aw. This indicates that in general rotation is sensitive to height 
change . 

In addition, the weight matrix M cannot be assumed to be dia- 
gonal but a full symmetric matrix as we shall see later on. Therefore, 
the two sets of constraints (2.4-11) and (2.4-12), in global spatial 
coordinates and (2.4-18) and (2.4-19) in local coordinates are a com- 
plete set as indicated in (2.4-7). This is the case regardless of which 
3-D measurement is being considered, and even if we are only interested 
in rotations. 

It is useful at this point to outline the approach taken by 
(Cannon, 1979; Cannon and Rochester, 1981) in establishing a terrestrial 
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reference frame by means of VLBI observatories since It leads to the same 
set of constraints. It Is based on analogies to the study of differen- 
tial deformations In continuum mechanics. 

Define the displacement (deformation) vector for baseline 1 

D = 3t - X (2.4-26) 

il 

This quantity Is In general not equal to zero due to deformations of 
the polyhedron and to measurement errors. An analog to the differen- 
tial tensor of the Infinitesimal displacement field Is given by the 
unitless matrix for the k'th baseline 


T 

1 0, 


(2.4-27) 


This matrix can be written as the sum of two matrices 


<=1 = i I «i 


where 


T 1 T T 

e^=C.+C. ^ (D,X^ +X^ DJ 


'1 ''1 ''1 


\|f 


(2.4-28) 


(2.4-29) 


is analogous to the strain tensor and describes deformations of the 
i'th baseline and 


T 1 'll >T’ 

= C -C. = ^ (D.X^ - X. D.) 

^ ^ ^ 'Ix^ f ' “i “i ^ 

1 


(2.4-30) 


is analogous to the rotation tensor of the displacement field. The 
displacement vector is then the sum of the "strain" matrix and the 
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"rotation" matrix times the initial polyhedron coordinate vector of 
station 1 


”1 = I ^ J \\ " 


( 2 . 4 - 31 ) 


A weighted mean global strain matrix is defined over the p 
polyhedron stations by 


E = 


P 

I m. 


P 

I me 
1=1 ^ 


i=l 


as well as a weighted mean global rotation matrix 


A = 


P 

Z 


Z m.n. 

1=1 ^ ^ 


i=l 


where m^ is a weighting factor. We can then define 


( 2 . 4 - 32 ) 


( 2 . 4 - 33 ) 


*1 = E + 


= A + (1)^ 


( 2 . 4 - 34 ) 

( 2 . 4 - 35 ) 


where and are residual deformations and rotations such that 


P 

Z m_. 


Z me 

i=l 


= 0 


( 2 . 4 - 36 ) 


i=l 
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1 


0 


(2.4-37) 


P 

Z m. 


P 

E m w . 
i=l ^ ^ 


i=l 


If the polyhedron stations were well distributed over the 
earth’s surface the matrix E would be representative of the global 
deformations of the earth. As can be easily shown, if trace (E) would 
differ significantly from zero (basically a weighted average of the 
distance changes) this would indicate global earth expansion or contrac- 
tion, l.e., a scale change. The off diagonal elements would indicate a 
global skewing ("shear") of the polyhedron. These effects could be 
absorbed by the CTS coordinates by computing new spanning base vectors 
as described by Cannon. The residuals £^, therefore should be due 
to non-global phenomena, l.e., deformations. 

Any global rotation of the stations would be indistinguishable 
from errors in the SNP transformation (i.e., in the estimation of polar 
motion and earth rotation or modeling of precession and nutation) . 
Therefore, the mean global rotation matrix should be identical to zero, 
i.e. , 

A = 0 (2.4-38) 

and thus (2.4-35) becomes 

= 0)^ (2.4-39) 

i.e., residual rotations should be almost entirely due to non-global 
phenomena . 
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Let us now examine the residual equations (2.4-36) and (2.4-37). 
They look quite similar to equations (2.4-11) and (2.4-12). If e^ and 
are tensor analogs, It should follow that so are and co^^. This can be 
seen from 



where 



0 

1 

0 


0 



(Y^ + Y^) 


1 



(2.4-40) 


(2.4-41) 


(2.4-42) 


Is the sum of a symmetric (e^) and antisymmetric (u^) matrix. The con- 
straints then become 


P 

Z m. 
1=1 


(Yi-YiHXi-X^)! 


0 


(2.4-43) 




(2.4-44) 


These equations are equivalent to (2.4-11) and (2.4-12) as well as to 
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(2.4-36) and (2.4-37). In the study of differential deformations in 
continuum mechanics, the differential motion of the particle displace- 
ment field can be split into two independent components, the strain ten- 
sor and the rotation tensor. It Is true that, in our case, we are 
considering small changes over baselines of possibly several thousand 
kilometer extent. However, the tensor analogy ends since as described 
above the 6 constraints form a complete set and are not independent. 
Furthermore, the M matrix is generally not diagonal. Therefore, our 
constraints should be written as 


^1 

^1 



T 

Y2 - Y2 • • • 
Y2 + Y2 • • • 


T 

Yp - Yp 
Yp + Yp 


M(X^ - Xq) 


0 


(2.4-45) 


or equivalently 


CMX = 0 


(2.4-7) 


Thus, we see that these constraints have appeared in the literature in 
different but usually less general forms (see also (Moritz, 1979; 

*r 

Richter, 1981)) particularly in the analysis of horizontal deformations 

i ^ 

(e.g., Brunner, et al, 1981). 


2.4.2 Reference Frame Constraints - Physical Approach 

In the previous section, we presented a set of constraints that 
insure that the reference frame defined through the fundamental poly- 
hedron is maintained, in an M-norm sense. In the following, we show 


42 


that this corresponds to maintaining a discrete analog of Tisserand's 
mean axes of crust. 

The time derivative of a position vector fixed in the earth is 
given by (e.g., Goldstein, 1981) 

Ml ■ Mt ^ 

The subscripts I and T indicate that the time derivatives are with 
respect to an inertial and terrestrial frame, respectively. The vector 
0 ) is the instantaneous angular velocity of the earth or its rotation 
vector. Equation (2.4-46) can be written as 


Vj. = + 0) X X (2.4-47) 

For the rigid body, = 0 since there is no rotation of points rela- 
tive to the terrestrial frame. In this case, results solely from 
the rotation of the earth in space. For the deformable earth 
denotes those motions relative to the chosen terrestrial frame, l.e., 
deformations. Furthermore, the rotation vectors are in general dif- 
ferent for each point, although for points on the same tectonic plate, 
they may be nearly the same. In any case, deviations in the rotation 
vector from one point to another are small and thus it is useful to 
define an instantaneous mean rotation vector, o) such that 

fff V • V pdE = minimum (2.4-48) 

g T T 

where p denotes density. This condition defines the motion of the 


43 


reference frame termed the Tisserand mean axes of body (Munk and 
MacDonald, 1960). According to (Smith, 1981), this integral should be 
evaluated only over the earth's crust (solid outer surface) since this 
is where our observations are taken and station coordinates assigned. 

In this case, one has a Tisserand *s mean axes of crust. Since observa- 
tions are only available at a finite number of observing stations (a 
polyhedron) condition (2.4-48) is unrealizable in practice. In this 
case, only a discrete analog of this condition is attainable 


P 

i (V * V ) = minimum • ' (2.4-49) 

1=1 i ^i 


where m^ are mass elements. Hopefully, the distribution and number of 
these stations will make this approximation meaningful, in the sense of 
being representative of the motion of the earth's surface. 

It is useful to present a set of constraints equivalent to 
(2.4-49). Without loss of generality we continue with summation instead 
of Integration for the reason mentioned above-. That is, we shall con- 
sider the polyhedron as a system of discrete mass particles with inter- 
nal motions and rotating in inertial space. Its angular momentum 
vector H is related to the torques L by Euler's equation 


L = 


dH 


dH 

dt 

I 

_dt_ 


+ u) X H 


(2.4-50) 


The total angular momentum is given by 


P 

H = Z m (X X V ) 
i=l 


(2.4-41) 
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which from (2.4-47) 


p 

H = E m [X X (w X X + V )] (2.4-52) 

i=l ^ i T 

P P 

= E m [X X (o) X X )] + E m (X x V ) (2.4-53) 

i=l ^ ^ i=l 

= I • 0) + h ' (2.4-54) 

Thus, the angular momentum splits into a rigid body motion 

= I • CO (2.4-55) 

where I is the Inertia tensor, and into a relative angular momentum, h 
due to the deformation of the system of particles. This is analogous to 
splitting the gravity field into a normal and disturbing potential, or 
a satellite orbit into a Keplerlan orbit and its perturbations. 

Carrying the satellite analogy further, the fundamental polyhedron cor- 
responds to the Keplerian orbit. If at the fundamental epoch, the 
earth became rigid, then the axes defined by the polyhedron would con- 
tinue to rotate with angular momentum H^. 

It can be shown that the constraints (2.4-48) are equivalent to 
the condition 

h «= 0 (2.4-56) 

Working along the lines of (Jeffreys, 1970) denote for the discrete case 
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(2.4-57) 


T = 


I m (V 

i=l 1 



P 


S m (V -0) X X )«(V -0) X X ) 
i=l ^ 

P 2 

2 [ (V^ - (ji^z + w^y) ^ 

+ - a)^x + 

+ (V^-oj^y + W2x)J] 


(2.4-58) 


(2.4-59) 


Minimizing T with respect to the components of to, the instantaneous 
rotation vector, yields 


3T 


P 2 2 

2 m, [ (y + z )w. 

i=l ^ 


+ (yv2-zv2) - x(o)2Z + W2y)]^ 


= 0 


(2.4-60) 


3T 

3(^2 


V 2 2 

Z m. [ (x + z )(jJ_ 
i=l ^ 


+ (xv^-zv^) - y(u)^x + u) 2 z)]^ 


= 0 


(2.4-61) 


3T 

3(03 


I m, [(x^ + y^)co + 
i=l ^ 


I- 

I- 

- z (tOj^x + ( 02 y) ] ^ = 0 


(2.4-62) 


In matrix form. 


P 

Z m. 
i=l ^ 


' 2 . 2 


“ 


“ “ 


— •“ 

y + z 

-xy 
2 , 2 

-xz 



P 

yv 3 - ZV 2 

-xy 

X + z 

-yz 
2 . 2 


<02 

= Z m. 
i=l ^ 

XV 3 - zv^ 

-xz 

-yz 

y + X 

i 

“3. 


_yv^ - xv^ 


(2.4-63) 
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or, 


I . 

implying that h = 0 from (2.4-54). This means that 



N 

1 

o 


Vn, 

p 



Tl 

= z 

2 0 -X 



i=l 



T2 


-y X 0 

i 

7 t 3 . 


N 

1 

o 


dx 

p 




= z 

N 

0 

1 


dy 

i=l 





-y X 0 

i 

dz 


(2.4-64) 


(2.4-65) 


(2.4-66) 


where dX^ are the differential displacements of the polyhedron stations 
with respect to the fundamental epoch. Considering that only periodic 
re-observations of the polyhedron deformations will be available over 
finite time intervals, the differential displacements can be approxi- 
mated by a finite displacement vector AX^. Finally, if we interpret the 
mass elements m^ as weights, the discrete approximation of the conditions 
(2.4-48) and (2.4-56) is equivalent to the set of orientation con- 
straints (2.4-11) of the previous section. As pointed out by (Munk and 
MacDonald, 1960) only the motion of the Tlsserand axes are defined by 
the above constraints, the choice of origin and orientation being arbi- 
trary. Following the requirements outlined in Chapter 1, we can choose 
the orientation of the polyhedron to be consistent with the BIH system 
and the origin to be at the center of mass, both at the fundamental 
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epoch. The translation constraints (2.4-12) maintain the origin 
definition. 

Thus, we have seen that the CTS reference frame as maintained 
by the constraints CMX = 0 Is a discrete Tlsserand's mean axes of crust 
or geographic axes as defined by (Munk and MacDonald, 1960). 

2.4.3 Alternate Approaches 

One drawback to Tlsserand axes Is precisely the problem of rela- 
tive motions as pointed out by (Moritz, 1979, 1980a). In the presence 
of secular motions the Tlsserand axes will rotate with respect to the 
observatories. In other words, the constraints CX = 0 (all stations 
are weighted equally, M = I) will Introduce Inconsistencies In the 
maintenance of the reference frame since In general the constraints have 
no relation to physical processes. This Is the situation If one Is 
restricted to geodetic observations on the earth's crust. 

It would be an Improvement If there was some model available 
for the expected deformations from which the model weight matrix M 
could be constructed. This can be done through absolute motion models 
from which plate motion velocities with respect to the mantle could be 
computed (section 4.2). Minster and Jordan (1978) refer to a mean meso- 
spheric frame which "Is fixed with respect to the average position of 
the deep mantle, assumed to be rigid or at least to have typical inter- 
nal motions with slower motion than the motion of the plates." They 
construct a model based on the Wllson-Morgan fixed hot spot hypothesis 
(Wilson, 1963, 1965; Morgan, 1971, 1972) which is used in the 


48 


simulations of Chapter 4. Another model is constructed by (Solomon and 
Sleep, 1974) by constraining the crust to have no net rotation with 
respect to the mantle. This is a requirement for any absolute motion 
model that is being considered for the CIS. It fits in with the 
generally accepted hypothesis that the crust and mantle very nearly 
rotate together (Smith, 1981). In any case, (Bender, 1981) indicates 
that the absolute plate motion velocities appear to differ by about 
1 cm/year among the available absolute motion models even though they 
are derived from different plate motion assumptions. It should be 
noted though that there is a controversy regarding the hot spot hypo- 
thesis and absolute reference frames (e.g., Le Plchon, 1973). 

Suppose a model weight matrix is constructed from an adopted 
absolute motion plate model (section 4.2). It can be constructed for 
any number of stations on the crust so one can define M_ as the global 
model matrix of infinite dimension. In this case 

CM_X = 0 (2.4-67) 

G 

for a model that does not contain any net rotations of the crust with 
respect to the mantle. Recall that X = are the coordinate 

changes (deformations) . For a particular polyhedron of stations a 
finite dimensional M is constructed, a subset of M^. Therefore, one 
could use the constraints CMX = 0 developed in the previous section as 
an alternative to CX = 0 when a model is available. 

Alternatively, one could use the constraints CX = Y where Y is 
computed from the model. If Y Q, this does not mean that there are 
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coimnon rotations (or translations) between the crust and mantle. 

Rather, for this particular distribution of stations, the absolute plate 
motion model Implies, for example, that the coordinate changes due to 
plate motions of the crust over the mantle do not sum to zero. Another 
approach would be to use the M matrix directly In the deformation 
analysis without Imposing any constraints. All these approaches have a 
corresponding deformation estimation algorithm which will be discussed 
In Chapter 3. Here we discuss In general terms the physical signifi- 
cance of using an M matrix. 

Recall the assumption that the mantle and crust rotate In a mean 
sense together, or equivalently that the M matrix Includes no net rota- 
tlons between the crust and mantle. The CTS reference frame axes, 
defined at the Initial epoch by the fundamental coordinates X^, can then 
be considered fixed with respect to the crust and mantle. At a later 
epoch, the CTS observatories will have different coordinates (where 
all global rotations and translations have been removed) with respect to 
this frame. These new coordinates In the crust and mantle fixed frame 
define the CTS. In order to maintain the reference frame consistently. 
It Is necessary to estimate the deformations X = X^ - X^. 

Summarizing, periodically all CTS stations observe In a short 
campaign for the purpose of monitoring deformations. From the baseline 
lengths and the model matrix M (If one Is adopted) , the deformation 
vector X Is estimated In a second adjustment using one of the four 
algorithms of the next chapter. Any Inconsistency between the 
re-estlmated baseline lengths and the plate model will show up here. 
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of course, the baseline lengths being crust fixed observations cannot be 
used to construct an absolute motion model. On the other hand, if there 
were no inconsistencies there would be no reason to challenge the 
assumption that the reference frame axes are crust and mantle fixed. 
After this short deformation analysis campaign, the subset of earth 
orientation monitoring stations continues its regular operations. But 
now, input into the earth orientation estimation algorithms (using the 
parameterization of section 2.3) are the updated CTS coordinates in the 
crust and mantle fixed system. Thus, polar motion and earth rotation 
still refer to the same reference frame axes as defined by X^, l.e., 
the ref erence system is maintained. If one did not correct for the 
deformations, the reference frame axes would rotate with respect to the 
CTS and degrade the earth orientation estimates. Applying the con- 
straints CX = 0 without an M matrix means that the reference frame is 
only fixed in the crust; and if the distribution of stations is inade- 
quate (CX 0 in reality) , then the earth orientation parameters would 
be degraded but to a lesser extent than not correcting for deformations 
at all. 

Of course, one could use an absolute motion model and not talk 
about the mantle at all. But implicit in M would be the mantle fixed 
frame and the assumption of no common rotations with the crust. In 
fact, it is possible to adopt only a relative plate motion model (e.g.. 
Minster and Jordan, 1974, 1978) and then no assumption about an absolute 
reference frame in the mantle would be involved at all. Everything 
would then refer to the crust, an approach that probably has many 
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advocates (e.g.. Smith, 1981). However, a relative motion model pro- 
vides information only, of course, on the relative motions between 
stations not on changes in station coordinates (deformations) . In any 
case, it seems preferable to estimate deformations using some plate 
model along with the re-estimated baseline lengths than to artificially 
impose the constraints CX = 0. This conjecture will be tested in the 
simulations of Chapter 4. It is important to estimate the deformations 
using the geophysical data in a weak way in line with, as stated in 
section 1.1, the requirement that the CTS should not be dependent on 
geophysical hypotheses. In other words, the deformation estimates 
should be as Insensitive as possible to errors in the geophysical model. 
Furthermore, any inconsistencies between the geodetic and geophysical 
data should be detectable in the estimation process. The purpose of 
the remainder of this investigation is to find such an estimation 
algorithm. 
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3 DEFORMATION ANALYSIS ALGORITHMS 


3.1 Introduction 

One of the requirements for the establishment of a reference 
system is the adoption of well-defined computational and estimation 
algorithms. We investigate in this chapter four possible algorithms 
for the analysis of polyhedron deformations. 

The analysis of deformation is accomplished in a two step pro- 
cedure. Periodically, observations are taken from all polyhedron sta- 
tions. Each measurement system analyzes its own data from which the 
estimated baseline lengths, along with their covariance matrices, are 
pooled into one common set. These baseline lengths are then compared 
to their corresponding values at the initial epoch which determines how 
the polyhedron has deformed. However, the absolute location of the 
deformed polyhedron, i.e. its new coordinates (or rather the change in 
coordinates relative to the initial epoch) , is undetermined from just 
the length of its edges, but this is what we seek. Without any other 
information the estimation of the defoxnnation vector in this form is 
singular due to the familiar origin and orientation defects and there- 
fore, the best linear unbiased estimate (BLUE) does not exist. This 
leads us to investigate other classes of estimators. Two are chosen 
from the class of biased estimators and two from the class of condi- 
tionally unbiased estimators. All are general enough to incorporate 


53 


adopted geophysical deformation models but do so in fundamentally dif- 
ferent ways. Furthermore, considering that a particular set of station 
coordinates defines the reference system, it is necessary to deal with 
the possibility of the loss or addition of a particular number of CTS 
stations. This leads to an application of least squares collocation. 

3.2 Mathematical Model 


Given are the adopted fundamental coordinates Xq and their cor- 
responding set of fundamental baseline lengths at an initial epoch. 
By comparing the estimated baseline lengths at a later epoch t to 
Lq, the deformation of the polyhedron can be estimated. 

The mathematical model for the deformation analysis is simply 
the chord length of baseline i - j at epoch t 

^*"ij " (3.2-1) 

It is fundamental that the linearization of this model be performed 
about the initial coordinates, X^, that define a datum for monitoring 
the time variations of the polyhedron. This is especially the case for 
our problem since we will be dealing with estimates that are sensitive 
to the initial parameter approximations. Linearization of (3.2-1) 
about Xq yields for baseline i - j 



(3.2-2) 
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Adding a true noise vector V yields the basic observation equation 
L = AX + V (3.2-3) 

The observation vector is 

L = - Lq (3.2-4) 

2 -1 

the changes in the baseline lengths. We denote a^P as the covari- 

2 

ance matrix of the observed baseline lengths where Oq is the variance 
of unit weight. The assumption is that any uncertainty in is not con- 
sidered (see section 4.4). The parameter vector X includes the 
deformations of the polyhedron, i.e., the change in coordinates between 

the initial epoch and a later one (e.g., AX , AX , ... of Fig. 1). 

^0 *^1 

The design matrix A contains the partial derivatives appearing in 
(3.2-2) evaluated at X^. It has dimensions n x u where n = p(p-l)/2 
is the number of polyhedron baselines (assuming all are observed) , 
u = 3p is the number of polyhedron coordinates, and p is the number 
of vertices. The column rank of A is deficient by 6, due to the 
orientation and origin defects so that 

R(A) = 3p - 6 (3.2-5) 

where R denotes the rank of the matrix. 

Without any a priori information on the parameter vector X 
one is limited to the Generalized Gauss-Markoff (GGM) estimation model 
(L, AX, Q^) where (e.g., Rao and Mitra, 1971) 

E{v} =0 ; D[V] = (3.2-6) 
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from which 


E{L} = AX ; D[L] = CqP"^ (3.2-7) 

The operator E stands for expectation and D for dispersion. 

The rank deficiency of the A matrix implies that the best 
linear unbiased estimate (BLUE) for X does not exist. This leads us 
initially to investigate the class of biased estimators. In the next 
section, we examine two estimators of this class. 

3.3 Biased Estimation 

3.3.1 Best Linear Minimum Bias Estimation 
3. 3. 1.1 A Deterministic Approach 
It is well known that the method of least squares can be 
developed in a purely deterministic manner, most easily using the con- 
cept of inner product spaces. The starting point is an inconsistent 
set of linear (in our case linearized) equations 

Y = AX (3.3-1) 

The matrix A represents a linear transformation of the vector X in 
e'^, a u-dimensional (parameter) vector space, to a vector Y in e'^, an 
n-dimensional (observation) vector space, e” becomes an inner product 
space with the definition of an inner product, in the most general case 
a weighted inner product 

<yi»y2>p = yi^y2 ’ ^l’^2^ (3.3-2) 
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It must fulfill the following properties (Davis, 1975) 


(1) 

<x,x>p ^ 0 , <x,x>p= 0 if and only is x = 

0 (Positivity) 

(2) 

<Xi,X2>p = <X2,x^>p 

(Symmetry) 

(3) 

< 0 Xj^,X 2 >p = a<x^,x^>p ; a real 

(Homogeneity) 

(4) 


(Linearity) 


For the weighted inner product to be properly defined, the weight 
matrix P must be positive definite. An inner product is essential 
for least squares solutions since it Introduces the concept of projec- 
tion. 

An inner product space is also a normed vector space, the 
weighted (ellipsoidal) norm being defined through the inner product as 

II y lip = ; yGE^ (3.3-3) 

providing the concept of length. It must fulfill the properties 
(Davis, 1975) 

(1) II y lip i 0 (Positivity) 

(2) II y lip ® 0 snd only if y = 0 (Definiteness) 

(3) II ay lip = Ia|||y||p for every scalar a (Homogeneity) 

(4) II x + y lip < ||x|| p+ ||y ||p (Triangle Inequality) 

Again, P must be positive definite. 

Normed spaces are also metric spaces, the weighted metric 
defined as 

d(yj^,y2)p “ II yi~y2 Up " (3.2-4) - 

providing a generalized concept of distance. 
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Likewise, in we define 


X 

= *1^2 ’ ^ 

(3.3-5) 

II X II “ (x'^^Mx)^'^^ ; X e e“ 

(3.3-6) 

d(Xi,X2)M = II x^-x^ llj^ = [(x^-X2)^M(Xj^-X2)]^'^^ 

(3.3-7) 


In the most general case, and the one that we will encounter 
in the analysis of deformations, the A matrix of equation (3.2-3) and 
the weight matrix M in (3.3-5) - (3.3-7) are rank deficient. For 
non-positive definite M and P , their corresponding weighted inner 
product and normed spaces are improperly defined as will be discussed 
later. In this investigation, though, P is always assumed to be 
positive definite. 

Consider a solution to the set of equations (3.3-1) 


X = GY ‘ (3.3-8) 

In our general case, we would like to preserve the property of least 
squares and this can be accomplished by minimizing the weighted norm 

|jY-AXl|p= [(Y-AX)^P(Y-AX)]^'^^ (3.3-9) 

This leads to the familiar normal equations 

NX - U = 0 (3.3-10) 




where 


N = 


T 

A PA 


U = 


T 

A PY 


(3.3-11) 
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Since A is in our case rank deficient, so is N . From (3.2-5) 

R(N) = 3p - 6 (3.3-12) 

and therefore (3.3-10) cannot be solved for X using Cayley inversion 
of N . In order to arrive at a solution of (3.3-10) , generalized 
matrix algebra is required, particularly that of the pseudoinverse so 
that the solution is unique, the primary prerequisite property. The 
ordinary pseudoinverse solution X = n"*”u where + denotes the pseudo- 
inverse is well known and has been used, for example, in the analysis of 
local horizontal deformations (e.g., Brunner et al., 1979). However, 
the most general case of a weighted M-norm, particularly for a singular 
M matrix, has not been treated let alone applied to any geodetic prob- 
lem. Let us, then, develop the case of an M-norm in e'^, considering 
first that M is positive definite, along the lines of (Rao and Mltra, 
1971). 

In order to Introduce the weighted norm |Ix]|„ we minimize the 

M 

Lagrangian function 

4) = X^MX - 2K^(NX-U) (3.3-13) 

so that the solution vector will have the additional property of mini- 
T 

mizlng X MX subject to the P least squares property implicit in the 
normal equations. Minimizing (J) with respect to X and the Lagrangian 
multiplier K yields two matrix equations 

MX - NK = 0 (3.3-U) 

NX - U = 0 (3.3-15) 
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From (3.3-14) 


X = (3.3-16) 

and substituting Into (3.3-15) 

NM“^NK - U = 0 (3.3-17) 

Since 

R(NM“^) = R(N) = 3p - 6 (3.3-18) 

then 

K = (NM“^N)‘^U (3.3-19) 

Substituting Into (3.3-16) yields the following unique solution 

X = m"^N(NM”^N)'*‘u (3.3-20) 

where for (3.3-8) 

G = m”^N(NM"^N)'^a'^P (3.3-21) 

Note that the + could be replaced by any generalized (g) Inverse 
(Appendix A.l). Using the symbolism of (Rao and Ultra, 1971) 

G - (3.3-22) 

is the minimum M-norm P least squares g-lnverse of A . It has also 
been referred to as the weighted pseudoinverse (Bouillon and Odell, 
1971). In the case of M = I, 


60 


(3.3-23) 


+ + T 

G = Api = N A P 

and 

X = n'^U (3.3-24) 

the ordinary pseudoinverse solution. 

In order to prove that G Is Indeed an A^^ the following four 
conditions must hold (Appendix A.l) 


AGA = A 

(3.3-25) 

GAG = G 

(3.3-26) 

(GA)'^M = MGA 

(3.3-27) 

(AG)**^? = PAG 

(3.3-28) 

Conditions (3.3-25) and (3.3-28) are equivalent to 

(see Appendix A.l) 

T T 

A PAG = A P 

(3.3-29) 

and (3.3-26) and (3.3-27) to 


T T 

G MGA = G M 

(3.3-30) 

The proof that (3.3-29) and (3.3-30) are fulfilled 
found In Appendix A. 2. 

It can be shown (Appendix A. 3) that 

for (3.3-21) can be 

= M"^N(NM“^)'^ 
IM 

(3.3-31) 


the corresponding G matrix for the consistent set of normal equations 


so that 
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(3.3-32) 




from which we get the relationship 


+ + T 

^>M = 


(3.3-33) 


Summarizing, the solution X = GY for the set of inconsistent 
linear equations (3.3-1) is unique and has the property of minimum 
M-norm P least squares. Substituting X = GY into (3.3-1) yields a new 
vector 


y' = AGY (3.3-34) 

where y’ = AX is a consistent set of equations. The matrix AG can be 
easily shown to be Idempotent and therefore a projection operator. It 
projects any vector Y in into y’ a vector in the image (the column 
space) of A . In fact, the consistency of Y* = AX is used in defining 
the concept of generalized inverse (Rao and Mitra, 1971) which indi- 
cates why it is important that be a g- inverse. This will concern 
us further in section 3. 3. 1.4 where the more general case of positive 
semidefinite M is treated. 

3. 3. 1.2 Estimation Model 

The results of the previous section can be applied to the 
linear estimation problem essentially through an appropriate choice of 
inner products by considering the weight matrices P and M as inverses 
of moment matrices. Consider again the set of observation equations 

L = AX + V (3.2-3) 


62 


We assume an expanded GGM model (L,AX,Q^,Q-) where 
E{V} = 0 

= E{w'^} = D[V] = CTqP"^ 

and 

Q- = E{3DC^} = (TqM ^ if Q- positive definite) 

- ^ 

such that 

y- = e{x> = X 

Z- = E{(X-X)(X-X)^} 

It follows, as in the GGM model (3. 2-6, 7), that 


(3.3-35) 

(3.3-36) 


(3.3-37) 


(3.3-38) 

(3.3-39) 


e{l} = AX ; D[L] ■= (3.3-40) 

In this setup, X is an Independent estimate of the parameter vector 

and Is stochastic In nature. On the other hand X is deterministic. 

In this section, we assume that Q- is positive definite. Note the two 

2 2 

variances of unit weight and T^, the latter related to the para- 
meter space. Note that X does not appear in (3.3-40), but enters only 
through the E^ inner product weight matrix M. 

The solution vector X (3.3-20) then becomes our first defor- 
mation estimate 
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= Q^(NQ^)‘^U (3.3-41) 

= m”^N(NM"^N)‘''u (3.3-42) 

= G^L = (3.3-43) 

GjL being defined as in (3.3-21). This is derived from minimization of 
(3.3-9) 

II Y-AX||^ . II V||^ - v'^PV (3.3-44) 

followed by, as before, 

II x|| = minimum (3.3-45) 

/V 


Thus, the statistical meaning of follows from the model defined by 

(3.3-35) - (3.3-40). The P and M matrices are defined now through 

the moment matrices and Qjj-. In our case, is the covariance matrix for 

the baseline lengths (since E{v} =0). In the following then, we 

denote a covariance matrix by Z and a moment matrix by Q , e.g., 

Instead of Q . The moment matrix Q- can be constructed from an adopted 
V A. 

model that predicts the deformations of the polyhedron stations (an 
example is given in section 4.2). This can be done only approximately 
for this estimation model since in (3.3-37) y- = X is unknown and we 
must approximate y- by X. 

Now that X^ has statistical meaning, using (3.3-40) 

E{X^} = Gj^E{L} = G^AX = M~^N(NM"^N)'^NX (3.3-46) 

But the rank of A is not full so that 
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(3.3-47) 


G^A 4 


I 

u u 


since in general for the product of two matrices A and B 


R(AB) £ R(A) and < R(B) (3.3-48) 

Thus, 

E(Xj^) X (3.3-49) 

and X^ is a biased estimate. However, it can be shown that X^ results 
from minimizing the bias norm 

||i-GA|| , (3.3-50a) 

or alternatively 

tr[(I-GA)M"^Cl-GA)'^] (3.3-50b) 

where tr denotes the trace operator (Chipman, 1964). For this to hold, 
it is necessary and sufficient that 


AGA = A 


(GA)^m"^ = m"^ga 


(3.3-51) 


An alternative formulation for the minimum bias estimator is as 

follows (Rao and Ultra, 1971). Considering the model (3.3-35)- 

T 

(3.3-40) find a linear estimate G L for X such that 

tr(SgTj^) = al tr (g\g) = Hg"^ 1|^ (3.3-52) 

Xj 

is a minimum (this is the minimum norm condition) in the class of 
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estimators which minimize the bias (the least squares condition) 


II aV - I |L = [tr((AV-I)\(AV-I))]^^^ 


(3.3-53) 


where Q- is positive definite. Then is a minimum (Z^^) -norm, 

T T 

-least squares solution of A G = I (which is inconsistent) such that 


- (a")* 


'’xV 
-1 -ll^ 


(3.3-5A) 


using a result from (Rao and Mltra, 1971) whereby 


and 




- G^L 


From the above, it can also be seen that is also a minimum variance 
estimate in the class of minimum bias estimators (see also (Chipman, 
1964)). 

The covariance matrix for Xj^ is given from (3.3-43) by 
,T 


= aQM"^N(NM~^N)'^N(NM”^N)'*'NM"^ 
2 + + T 


(3.3-55) 


where is given by (3.3-31). An unbiased estimate for the a posteriori 
variance of unit weight is (Uotila, 1967) 


66 


(3.3-56) 


..X - 

V PV 


where 


-2 = 

°0 n - R(A) 


^ 'P aT’ 

V PV = L PL - xju 


(3.3-57) 


In our case, for p stations 

„ .£ipLI 

and 

R(A) = 3p - 6 
from which it follows that 


n _ R(A) = (p-4) (p -31 = 1 ( 2 .- 3) !. 

" 2 2 (p-5)! 


(3.3-58) 


(3.3-59) 


(3.3-60) 


It is evident that in order to achieve redundancy at least five poly- 
hedron stations are required. 

2 

An estimate for Tq of (3.3-37) is approximately given by 


aT a 

^2 _ X MX 
”^0 u 


(3.3-61) 


where, in this case, u is the number of parameters (and the rank of M) . 

Summarizing, X^^ has the deterministic properties of minimum 
M-norm in the class of P least squares estimators and the stochastic 
properties of minimum variance in the class of minimum bias estimators. 
Therefore, Xj^ is known as the Best (minimum variance) Linear Minimum 
Bias Estimate or BLIMBE (Rao and Mitra, 1971). In the next section, 

A 

we see that Xj^ can be computed by augmenting the singular normal equations 
(3.3-10) with the set of constraints CMX = 0 derived in Section 2.4.1. 
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3. 3. 1.3 The Weighted Inner Constraint Estimate 
It is well known that the ordinary pseudoinverse estimate 



(3.3-62) 


is equivalent to augmenting the system of singular normal equations 
NX = U by a matrix C such that the conditions 


T 

NC = 0 


(or AC = 0) 


(3.3-63) 


CX = 0 


(3.3-64a) 


are fulfilled (Melssl, 1969; Blaha, 1971). The set (3.3-64a) is known 
as an inner constraint. For our particular case, we need six such con- 
straints, the number of rank deficiencies of N . We generalize the 
results to deal with the definition of an M-norm in the parameter 
space. We assume that M is positive definite. The derivation of 
this weighted inner constraint estimate is a generalization of the 
development in (Blaha, 1971). Instead of (3.3-64a) we substitute the 
condition 

CMX = 0 (3.3-64b) 

We start by minimizing the Lagranglan function 

4> = v'^PV - 2K^(AX + V-L) -2K2(CMX) (3.3-65) 

with respect to the unknowns V , X , , ^2 ’ some matrix 

manipulations we arrive at 
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N MC 


CM 0 -K 


(3.3-66) 


The resulting augmented normal matrix Is non-singular and we could 
solve for X directly. In order to eliminate K 2 > and derive an 
explicit expression for X we can write 


N MC I>2 


CM 0 D-D, 
3 4 


I 0 


0 I 


Explicitly, this represents four matrix equations 


+ MC = I 


ND- + MC D, * 0 
2 4 


CMDj^ « 0 


CMD 2 = I 


Pre-mult Ip lying (3.3-69) by C 


CND- + CMC D, = 0 
2 4 


But CN = 0 from (3.3-63) and since CMC Is full rank (which follows 
from the assumption of M positive definite) , we have 


^ = 0 


Pre-multlplylng (3.3-68) by C 


CND^ + CMC = C 
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from which 


= (CMC^)"^C (3.3-75) 

Inserting D 2 = into (3.3-71) 

(CMC^) (CMC^)"^ = I 
which implies that 

©2 = r>2 = C(CMC'*’)“^ (3.3-76) 

Substituting into (3.3-68) yields 

ND^ + MC'^(CMC^)”^C = I (3.3-77) 

This is equivalent to 

[N+kMC^CM][D^+k"V(CMC^)"^(CMC^)"^C] = I (3.3-78) 

This follows from (3.3-63) and (3.3-70). The scale factor k is 

arbitrary as can be seen by performing the multiplications on the left 

hand side of (3.3-78). It is useful for improving the condition of 
T 

N + MC CM for inversion. It can be shown following the same reasoning 

T 

as in (Blaha, 1971) that the matrix N + kMC CM is non-singular. There- 
fore 

= (N+kMC^CM)"^-k~^C^(CMc’^)"^(CMC^)“^C (3.3-79) 

From (3.3-66) we have 


70 


(3.3-80) 


X = D^U 


= [ (N+kMC^CM)"^ - k"^c'^(CMC^)“^(CMc'^)‘^C]U 


= (N+kMc’’‘^CM)"^0 


(3.3-81) 

(3.3-82) 


since CA =0, and 


E- = aQ[(N+kMc'^CM)"^ - k ^C’^(CMC'^)~^(CMC’^)"^C] 


(3.3-83) 


In order to demonstrate that this estimate is equivalent to 

+ T 

the BLIMBE estimate = Ap^^L it is necessary to prove that G = D^A P 

is For this to hold, it must satisfy the conditions (3.3-29) 

and (3.3-30). First, we show that 


T T 

A PAG = A P 


(3.3-29) 


Post-multiply (3.3-77) by A P 


T T T T —1 T T 

A PAD^A P + MC (CMC ) CA P = A P 


(3.3-84) 


which follows from CA =0. Second, we prove that 


T T 

G MGA = G M 


(3.3-30) 


Transpose (3.3-77) and pre-multiply by (D^a'^P)'^M (=G^M) 


TT T TTT T— 1 

(Dj^A P) MD^A PA + (Dj^A P) MC (CMC ) CM 

= G^MGA + PADj^MC^(CMC^)”^CM 


(3.3-85) 


T 

= G M 
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which follows from (3.3-70). Therefore, the weighted inner constraint 
estimate (3.3-82) is equivalent to the BLIMBE (3.3-43). However, the 
former is computationally more efficient considering that only a 
Cayley Inverse is required. 

A note should be added concerning the units for the computation 
of X by (3.3-82). In order for (3.3-64b) to be consistent, C and 
M must be unitless so that the 0-vector will be in units of length. 
Therefore, we divide the length unit elements of C by a mean earth 
radius R , i.e. (see 3.3-90 below). 



In addition, M is divided by its Euclidean norm 

II mH = (Zm^j)^^^ C3.3-87) 

so that 


5 = lifii (3.3-88) 

Therefore, in (3.3-82) we should write strictly 

X = (N+MC^(M)"^U (3.3-89) 

but this has no effect on the estimate which is independent of a scale 
factor k as shown above. 
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The constraints C 


from the condition AC = 0. 


C = 


••• 


I I ... I 


for this particular problem can be derived 
It is found that 

(3.3-90) 


fulfills this condition which is just (2.4-9), the reference frame 
maintenance constraints. Therefore, the BLIMBE (or equivalently the 
weighted inner constraint) estimate provides an algorithm to maintain 
a discrete version of the Tisserand mean axes of crust which incorporates 
also the stochastic properties of non-perfect measurements and adopted 
geophysical models (that come through the weighted parameter norm) . In 
this case, we can say that the reference frame axes are maintained in a 
minimum M-norm P least squares and BLIMBE sense by a specified number 
of CTS stations. 


3. 3. 1.4 A Generalization for M-Seminorms 
In the previous sections j we have derived X under the assump- 
tion that P and M are positive definite, i.e.. 


II X 11^ = (x'^Mx)^^^ > 0 for all xGe“ (3.3-91) 

II y lip = > o for all yGE^ (3.3-92) 

except for the trivial cases x = 0 and y m o for which equality holds. 

As we shall see, it is quite possible in our applications that the 
moment matrix Q- is only positive semideflnite (it must be at least 
positive semi-definite since it is a moment matrix). This more general 
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case raises several problems. First, M can no longer be defined as 
the Cayley Inverse of Os general, non-posltlve definiteness does 

A 

not necessarily Imply singularly, however, since we are dealing with 
moment matrices, non-posltlve definiteness does Imply singularity). 
Second, assuming an appropriate M could be found, the properties of 
the weighted Inner product and norm listed In section 3. 3. 1.1 may no 
longer hold for 

If M Is a positive definite matrix, It can be expressed as 

M = UAu"^ C3.3-93) 

by singular value decomposition where A Is a diagonal matrix of dimen- 
sion u whose non-zero elements are the real eigenvalues of M , 

T 

and the columns of the n-dlmenslonal orthogonal matrix U (U U = I, 

T 

UU = I) are the corresponding normalized eigenvectors U^. In this 
case, all the eigenvalues are positive (nonzero). The eigenvectors 
form an orthogonal basis for e'*. Therefore, the weighted norm 
II X 11^ and inner product <x^,x| can be defined properly over the 

entire space e'^. However, if M is only positive semidef inlte, then 
it can be decomposed as (e.g., Lanczos, 1961) 

M = U A (3.3-94) 

UU uppppu 

where A is a p-dimensional diagonal matrix whose diagonal elements are 
P 

the p non-zero eigenvalues of M (there are u-p zero eigen- 
values) and the p columns of the semiorthogonal matrix U 
~T~ ~~X 

(U U = I, UU I) which contain their corresponding p eigenvectors. 
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l.e., the reduced elgenspace 

[(X^.U^) I > 0 , 1 = (3.3-95) 

The are also called the principal vectors of M (Ben Israel and 
Greville, 1974). In this case, the principal vectors form on ortho- 
gonal basis only for a subspace of e'^, call it e'^. 

It follows that the weighted (ellipsoidal) norm and inner 
product can be defined properly over the same subspace E . That is, 
for the inner product 

<Xj^,X2>^ = x][mx2 ; e“ (3.3-96) 

and for the norm 

II X » (x'^Mx)^^^ ; xE e" (3.3-97) 

However, the parameter vector X may not necessarily be an element of 
e'^. Therefore, one is led to define an "improper" weighted (hyper- 
bolic) norm for the entire space E*^ although it can "properly" be 
applied to physical problems (Pease, 1965). Property 1 of a proper 
inner product (see section 3. 3. 1.1) must be modified to 

(l’) <x,x>„ > 0 (Positive Semidefiniteness) 

since the inner product could be equal to zero for x 0. This means 
that a vector x could conceivably be orthogonal to Itself but in the 
semideflnlte weighted sense. Likewise, the property (2) of definite- 
ness no longer holds for a weighted norm. Rao and Mitra (1971) refer 
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to such norms as seminorms. Furthermore, as shown in (Pease, 1975), 
the triangle inequality property (4) of section 3. 3. 1.1 may also fail. 
Consider as only positive semidef inite. Then 

as described above. In this case, a reasonable choice for M is 

M = (3.3-99) 


a well known result for the computation of a pseudoinverse. It can be 
shown that M is also positive semidefinlte (Lewis and Newman, 1968). 


Consider Q- as an operator 

A 

Ui U2 

Qx 5 E E ^ 

""l “2 

where^E and E are p-dimensional inner product spaces. For 
u ^ E ^ 


(3.3-100) 


V = ° 

Similarly, 

M : E E 


(3.3-101) 


(3.3-102) 


This holds since M = Q is a unique operator which makes the above 
choice of M suitable for our purposes (in fact, a reflexive g- inverse 
(Appendix A.l) also has this property). In general though, for any 
matrix A 
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which makes just a g- inverse unsuitable. 

Now that we have defined a seminorm 1| X 1|^ and found a suitable 
choice for M, let us address the problem of generalizing the deforma- 
tion estimate (3.3-43). This discussion will be primarily based on 
(Rao and Mltra, 1971; Mltra and Rao, 1974). 

Let us first consider the simplest case of an ordinary least 
squares estimate (i.e. N is full rank) with a weighted parameter norm 
defined in the parameter space where M is positive definite. We 
start from the normal equations 

NX - U = 0 (3.3-15) 

and minimize 

4) = X^MX - 2K^(NX-U) 

with respect to X and K . This yields the relationship 
MX - NK = 0 

in addition to (3.3-15). It follows that 

X = m"^ NK (3.3-105) 

Into (3.3-10) and solving for K 

K = (NM“^N)“^U (3.3-106) 


(3.3-103) 


(3.3-104) 


Substituting into (3.3-105) gives 


(3.3-107) 


X = m"^(nm"^)"^u 


= n“^u 


Of course, this could have been directly obtained from (3.3-5). How- 
ever, it does prove that the ordinary least squares estimate is invari- 
ant with respect to any defined positive definite norm in E^. This is 
not surprising since both N and M are full rank and are not 
restricted to any subspace of This is really what we mean by an 

A. A IX 

unbiased estimate since E(X) = X can conceivably hold for any XGE . 

In fact, it is enough that N be full rank even if M is not positive 
definite. That is, the best linear unbiased estimate (if it exists) is 
unaffected by any weighted norm in the parameter space. This follows 
from the fact that I - GA = 0 so that the bias norm (3.3-50) is zero 
independent of any weight matrix. 

A more general case will now be described where the M matrix 
is positive semldefinlte and the design matrix A is rank deficient. 
As pointed out in (Rao and Mitra, 1971) care must now be taken since 
both M and N are not of full rank. Let us return to the two matrix 
equations 


MX - NK = 0 


NX - U = 0 


(3.3-14) 

(3.3-15) 


Combining (3.3-14) and (3.3-15) 


(N + M)X - NK = U 


(3.3-109) 
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from which 


X = (N+M)"^(NK+U) (3.3-110) 

assuming that N4M is full rank the necessity of which for our applica- 
tion will be described below. Into (3.3-15) yields 

N(N+M)"^(NK+U) = U (3.3-111) 

that can be solved for K as 

K = [N(N+M)"^N]‘‘‘[I-N(N+M)“^]U (3.3-112) 

Substituting back into (3.3-110) 

X = (N+M)"^[N(N4M)“^]‘^U (3.3-113) 

- (N41tf)"^N[N(lHM)“^]‘’'N(N+M)“^U 
+ (N4M)“^U 

The last two terms can be shown to cancel using the relations (Rao and 
Mltra, 1971) 

A(A^PA)®A pa = a (3.3-114) 

(a’^PA)(a'^PA)®A^ = A*^ (3.3-115) 

T 

for any matrix P such that R(A PA) = R(A). This holds automatically 
for positive definite P for which the above results could be modified 
as 
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PA(a'^PA)®A^PA = PA 


(3.3-116) 


(A^PA)(A^PA)®a''^P = A^P (3.3-117) 

Using (3.3-117) the second term in (3.3-113) can be written as 
(N4M) [N (N4W) (N4M) "^NN^A^^PL 

which reduces from (3.3-114) to 
(N+M)~^NN®A^PL 


Applying (3.3-117) again, the second term is seen to be equal to the 
third term in (3.3-113) so that 


X = (N4M)"^N[N(N-«I)''^]‘‘'U 
= GL 


(3.3-118) 


where 


G = (N-W0"^N[N(N+M) ^N]‘*’a'^P 


(3.3-119a) 


Note that the same estimation model is assumed as in section 3. 3. 1.2. 
The relation (3.3-117) can be used to show that X is invariant with 
respect to any g- inverse for the term in brackets. It gives 

G = (N4M)"^N[N(N+M)“^N]®NN^a'^P (3. 3-1 19b) 

Another result from (Rao and Mitra, 1971) is that A(a"^PA)®A^ Is invari- 
ant for any choice of (A*^PA)® which applied to N[N(N+M) proves 

the above assertion. 
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According to (Mltra and Rao, 1974) for the possibly inconsistent 
set of linear equations L = AX 

G = (N+M)®N[N(N-fW)®N]®A^P 

called an Ap^ inverse of A, is one choice for a minimum M-seminorm 
P-semileast squares inverse of A . It follows that X = GL has minimum 
M-seminorm among the semileast squares solutions of L = AX. According 
to the same theorem, Ap^^ is unique if and only if N + M is positive 
definite explaining our use of Cayley inversion in the above derivation 
for X . However, uniqueness is not sufficient for our purposes since 
Apj^ may not even be a g-inverse (AA^^ A) . We require this property 
in order to make (3.3-34) 

l' = AGL 


consistent. Also, we would like A^^ to be ref lexive so that 
R(Apj^) = R(A) (Rao and Mitra, 1971). If these properties hold, then 
Apjj is called an but recall that unlike the situation in section 
3. 3. 1.1, P and M may be positive semldeflnlte. If both are positive 


definite, A^^^ is equivalent to the G matrix of 


BLIMBE. 


For Apjj to exist it is necessary and sufficient that 


C(M) n C(a'^) C 


(3.3-120a) 


where C denotes column space, according to theorem 3.6 of (Mitra and 
Rao, 1974). In our application P is always positive definite so 
that this condition could be modified to 
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C(M) n CCA*^) c 


(3.3-120b) 


which always holds. Therefore, in order for G to be Ap^, according 
to theorem 3.2 in the same paper, the following four conditions are 
necessary and sufficient 


PAGA = PA 

(3.3-121) 

MGAG = MG 

(3.3-122) 

(GA)*^M = MGA 

(3.3-123) 

(AG)"^? = PAG 

(3.3-124) 


That G fulfills these four conditions is proven in Appendix A. 4. 

For P positive definite (3.3-121) is equivalent to (3.3-25). Condi- 
tions (3.3-123) and (3.3-124) are equivalent to (3.3-27) and (3.3-28), 
respectively. Therefore, the difference between the M positive 
definite case of section 3. 3. 1.2 and M positive semidefinlte is just 
condition (3.3-122). Unfortunately, without the fulfillment of (3.3-26) 

A 

which can occur only if M is positive definite, the solution X = GL 
for G in (3.3-119) although M-seminorm P least squares is no longer 
BLIMBE as we shall see in the simulations of Chapter 4. Thus, we will 
call = Gj^L a MINDLESS (minimxim norm (seminorm in this case) least 

squares solution) using the terminology of (Schaffrin, 1975)* For posl- 

_ 

tive definite M , X^ = X^ as shown in 3.3. 1.2. Nevertheless, the 
solution (3.3-119) may still be applicable to deformation analysis as 
will be tested in Chapter 4. 
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The covariance matrix of 


^ is, using (3.3-118), 


Z. 

X 


where Oq 


* = GZ 
1 

= Oq (N+M) " [N (N-rtl) ”^N ] "^N [N (N+M) ”^N ] "^N (N4M) 


(3.3-126) 


is given by (3.3-56). 


Here, Tq is given approximately by 


;2 _ X^MX 
^0 “ R(M) 


where = Tq (compare to (3.3-37) and (3.3-61)). 


(3.3-127) 


The deformation estimate X^ of this section is no longer equi- 
valent to the weighted inner constraint solution of section 3. 2. 1.3 
(recall that there M was assumed positive definite). Therefore the 
most one can say is that the reference system is maintained in a mini- 
mum M-seminorm P-least squares (MINDLESS not BLIMBE) sense by a speci- 
fied set of CTS stations. 


3.3.2 Best Linear Estimation 


Another possible biased estimator is called the Best Linear 
Estimate or the BLE although several versions are available in the 
literature (e.g. Rao, 1973). Recall the multivariate definition of the 
mean square error 


MSE(X) = e{(X-X)(X-X)'''^} 

= Z- + [X-E{X}][X-E{X}]'^ 


(3.3-128) 


as a sum of covariance and bias squared. As shown in section 3.3.1, the 
BLIMBE minimizes the bias term. Conditional on this property is the 
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minimum variance one. The BLE, on the other hand, minimizes the mean 
square error. Unlike the BLIMBE, It can only be approached In a prob- 
abilistic manner. Its existence Is dependent on some a priori knowl- 
edge of the parameter vector as we shall see below. 

As before we will derive a homogeneous estimate of the form 

X = GL 

We assume now that both X and L are random variables with moment 
matrices 

= E{XX'^} (3.3-129) 

= E{Ll'“^} (3.3-130) 

Q^=»E{Xl’^> (3.3-131) 

where Q_ Is non-singular (note that this Is not assumed for Q ) . By 
the Gauss-Markoff theorem, the linear minimum mean square error esti- 
mate of X Is (Llebelt, 1967) 

^ (3.3-132) 

with mean square error matrix 

- <’x - (3.3-133) 

Note that some knowledge of the moment matrices Is required. 

For the observation equations 

L = AX + V C3.2-3) 
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we assxime an estimation model of the form (L,AX,Q^,Q^) where X and V 
have a priori probability distributions. For "V 


e{v} = 0 (3.3-134) 

= E{w'^} = D[V] = (3.3-135) 


and for X 

e{x} = X (3.3-136) 

D[X] = E{(X- X)(X-X)'^}= (3.3-137) 

(recalling that X is an independent estimate of X) so that 

= E{xx'^} = Zjj + xx"^ = Oq (3.3-138) 

2 + 

where for positive semidefinite Q^, M . Here X is stochastic 

and X deterministic, the oppositve of the BLIMBE model. The conditional 
distributions of L are then given by (Chipman, 1964) 


E{l|V} = AX + V ; D[L|V] = 

E{L|X} = AX ; D[L|X] = 
from which 

e{l> = e{e{l|v}> = e{ax+v> = AX 

D[L] = E{D[L[V]} + D[E{l|v}] 

= e{AZ^'^} + D[AX + V] 

T 2-1 

- V - =L 

Therefore, 


(3.3-139) 

(3.3-140) 


(3.3-141) 


(3.3-142) 
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= E{LL^} = E{ (AX + V) (AX + V)^} 

= A A^ + Oq p"^ (3.3-143) 

( T T 2-1 

= A A + A X X + Oq P 

= \+ E{L} E{L^} 

Furthermore, 

= E{XL^} = E{X(AX + V)^} 

= Eixx"*^} A^ + E{XV^} 

= (3.3-144) 

where we have assumed for both (3.3-143) and (3.3-144) that 

= E{XV^} = 0 (3.3-145) 

In our context, this means that baseline measurement errors and deforma- 
tions are uncorrelated. As we shall see, it is useful to assume as 
indicated in (3.3-138) that Q.. is also known to within a scale factor 

A 

2 

Oq, which can be taken as unity. The deviation of an unbiased estimate 
2 

for Oq, derived later, from unity will indicate the degree of compatibility 
between the baseline measurements and a chosen geophysical deformation 
model. With this in mind and from (3.3-132) it follows that 

Xg = Qx L (3.3-146) 

where S denotes that X is assumed stochastic. The mean square error matrix 
can be computed from (3.3-133) as 

MSE(Xg) = Oq [Q^ - + P“S"^ A (3.3-147) 

Note that strictly Q should be replaced by in (3.3-146) and (3.3-147) 

A 

although we will continue with this notation since using m"*" may indicate 
misleadingly that pseudoinversion is required. 
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We see that this estimate is general enough to deal with rank deficient 

A and Q matrices encountered in the analysis of deformations. In the 

case that Q is nonsingular, using the identities (Liebelt, 1967) 

X 

= (a'^PA+Q^"^)"^a'^P (3.3-148) 

and 

(a'^PA+Q^"^)"^ = Qx" P“^)"V (3.3-149) 

both for possibly rank deficient A , we get in simpler form 

X = (N+M)“^U (3.3-150) 

s 

and 

MSE(Xg) =aQ(N+M)"^ (3.3-151) 

where M = Q 

X 

We see that the approach here to overcoming the singularity 

problem due to the rank deficiency of A (and N ) is to use a priori 

information for the deformation vector X . In BLIMBE, the approach 
was to use pseudoinverse algebra or equivalently to augment N with a 
set of inner constraints. 

Now, it is quite useful to derive another estimate which is a 
limiting case of the above estimate Xg, following essentially (Bibby 
and Toutenberg, 1977). We minimize the quadratic loss function 

R(X) = e{(X-X)^B(X-X)} (3.3-152) 

The expression in the brackets represents a generalized distance (com- 
pare to (3.3-7) for the metric defined by B (assuming that B is 
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non-singular) , or a weighted norm 

II e II , ; e = X - X 

B 

where B indicates the relative importance of the 

2 -1 

of X . We start with the GGM model (L, AX, (T^P 
L = AX + V ; V ~ (O.OqP"^) 

The estimate X will be a function of L 
X = GL = GAX + GV 
from which 

X - X = (GA- I)X + GV 


and 


R(X) = E{[(GA-I)X + GV]’’^B[(GA- DX+GV]} 
= x’^CA^^- I)B(GA-I)X + e{v'^GBg'^V} 


where it can be shown that 


e{v'^GBG^V> = Og trace (BGP~^G^) 

2 -1 T '' 

Note that cTqGP G is the covariance matrix of X . 

with respect to G yields 


T T T T 2-1-1 
G = XX A (AXX A + O^P 


(3.3-153) 

different elements 
of section 3.2 

(3.3-154) 

(3.3-155) 

(3.3-156) 


(3.3-157) 


(3.3-158) 
Minimizing R(X) 

(3.3-159) 
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so that 


T T T T 2 —1 —1 

Xjj = XX A (AXX A +OqP ) L (3.3-160) 

where D denotes that X Is assumed deterministic. Note that this 
result is independent of B . The estimate doesn't seem to be of any 
use since it depends upon the true value of X which, of course, we do 
not know. 

Let us examine (3.3-160) more closely. By the identity 
(3.3-149) 

[AXx'^^A'^ + aQP"^]"^ = [P-PAX(x’^a’’^PAX+0q)"Va'^P] (3.3-161) 

°0 

Note that although this identity holds also for a rank deficient A 
matrix, in the literature full rank is usually assumed in the deriva- 
tions of BLE’s. Substituting into (3.3-160) 




Xx\'^[P- PAX(x'^NX+aQ)"^x'^A'^P]L 

'"o 


5" 


1 - 


T 

X^NX 


T 2 
X NX+Oq 


T 

XU 


X 


T 

XU 


x"^NX+a 


2 

0 


Now, 


(3.3-162) 
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X 


E{iL} = -1~ 2 X"E(U) 

X NX+a„ 


x'^NX 
x'^NX + ^0 


(3.3-163) 


2 1-1 


= X 1 + 


T 

X NX 


which shows that on the average,^ (3.3-160) Is an underestimate of X . 

Let us then compare the estimates (3.3-146) and (3.3-160). 

First, Xg is derived under the assumption that X is stochastic while 

A, ^ T 

for Xj^ it is assumed deterministic. Thus, for X^, XX is required 

/V r Ti 

while for X only its expected value, i.e. the moment matrix Q = E{XX }. 

Furthermore, X^ is on the average an underestimate of X. 

How can one alleviate the unfortunate situation of X^ for which 

the true value X is required. Several approaches are possible. If 

there is a priori information on X , say X, one could replace X by X in 

(3.3-160). This is essentially the approach advocated by (Rao, 1973). 

In practice, a covariance matrix may be available for X . There- 

T 

fore, a preferable approach is to replace XX in (3.3-160) by its expec- 
tation Q = e{XX } given in (3.3-138). This means using X instead of 

A D' 

Xj^. Therefore, we shall refer to X^^ as our Best Linear Estimate (BLE) . 

Consider the limiting case when the expected deformations X 
are known perfectly. Then Z— = 0, X = X and X^ becomes X^ which 
yields, as shown above, an underestimate of X . On the other 
extreme if Z- -► “, it is as if one has no a priori information on X so 

Jv 

/s 

that neither Xg nor X^^ are useful. In this case, only BLIMBE is 
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2 2 
available with M = k I (it is invariant with respect to k as is easily 

seen by examining (3.3-42)). We can then write the I-norm BLIMBE or the 

ordinary pseudoinverse solution as 


X = n'^U = lim [(N + k^I)~^U] 
k^-K) 


(3.3-164) 


for singular N. Then k = 0 In the case of no a priori information. Or 
more generally it is a limit for BLE when M 0. The expression in 
brackets in (3.3-164) is a special case of BLE, called the ridge esti- 
mator (see (Pavlis, 1979) for a good review). 

Considering these two limiting cases it should follow that 


^LE ^LIMBE^I 


(3.3-165) 


as we shall show below. Recall that the minimum norm property of 
BLIMBE is conditional on P-least squares. 

Consider the minimization of the Lagranglan function 

(j) = v"^PV + x'^MX - 2K‘(AX + V-L) (3.3-166) 

with respect to X , V , K which gives (we assume here that M is 
positive definite) 


PV - K = 0 (3.3-167) 

MX - a\ = 0 (3.3-168) 

AX + V - L = 0 (3.3-169) 

Solving (3.3-167) for K and substituting into (3.3-168) gives 
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MX - a'^PV = 0 


C3. 3-170) 


Substituting V from (3.3-169) into (3.3-170) yields 
X = (N + M)”^U 

as in (3.3-150). Therefore, the BLE has the property of 
V^PV + X^'^MX = minimum 

Denote 

X^ = BLIMBE 
X 2 = BLE 

As shown by (Hoerl and Kennard, 1970) 

(L - AX2)^ P(L - AX 2 ) 

= (L - AX^)^ P(L - AX^) + H 
H = (X2-Xj^)^N(X2-X^> 

Since H is a quadratic form and in general X^ f X 2 (see (3 
below) then 

(V*^PV)2 - (v’^PV)j^ > 0 

which already follows from the P-least squares property of X. 
(3.3-166) 

(v'^PV)2 + (x’^MX) 2 < (V*^PV)j^ + (X*^MX)^ 


(3.3-171) 


(3.3-172) 


3-175) 


(3.3-173) 
. From 


(3.3-174a) 
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or 

0 < (v'*^PV)2 - (v'^PV)^ < (ic’^MX)j^ - (x'^MX) 2 (3.3-174b) 

which proves (3.3-165), i.e. , between BLE and BLIMBE, BLE has the mini- 
mum norm property. Therefore, the BLE for the deformation vector is on 
the average smaller than the BLIMBE, it is closer to zero. On the 
other hand, the BLIMBE has minimiam bias. Mote that this comparison has 
been performed for only M positive definite recalling the discussion 
of section 3. 3. 1.4 for M positive semideflnite. 

A, A 

Since the vector U is the same for and X 2 it easily fol- 
lows that 


X 2 = (N+M)"^X^ (3.3-175) 

A A 

indicating that in general X^^ f X 2 , in the case when N is singular, 
except as explained before if M=0 in which case 


A + ^ 

X 2 = N NX^ = N U = X^ 


As a summary, it can be shown (Chlpman, 1964; Rao, 1976) that 
the BLE can be obtained from minimizing the mean square error matrix 
(3.3-128) also called the risk function 


R(G) = (I- GA)Q^(I-GA)'^ + 

A L 

2 

= (bias) + covariance 


(3.3-176) 


This yields 
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(3.3-146) 


X 2 = (N4M) ^U (if Q^ non-singular) 
as before. The minimum risk then becomes 
R(G2) = (I-G2A)Q^ 


(3.3-150) 


(3.3-177) 


which is just the mean square error matrix (3.3-133). 

It is useful to compute an expression for V PV. In general. 


<^T ^ ^ T ^ 

V PV s (L- AX) P(L- AX) 

T ^ 

*= L PL - 2X U + X NX 

rn ^ ^ 

= L PL - X U + (X N- U )X 


(3.3-178) 

(3.3-179) 

(3.3-180) 


V/hen Q is positive definite and from (3.3—150) 

A 


or 


and 


(N-rtI)X2 - U = 0 

A»n A»r 

X 2 N + X 2 M - U = 0 

aT T ^ ^ 

(X2N-U^)X2 = -^ 2^2 


Substituting into (3.3.180) yields 

aT T /vT ^ 

V 2 PV 2 = L PL - X 2 U - X 2 MX 2 


(3.3-181) 


(3.3-182) 

(3.3-183) 


(3.3-184) 
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From ( 3 . 3 - 166 ) it follows that the BLE minimizes (Q positive definite) 


V2PV2 + X^MX^ = L PL - U 


( 3 . 3 - 185 ) 


Now, 


E{V 2 P ^2 + X2MX2} = E{L^PL} - E{X 2 U} 


Since 


E{l'^PL} = E{tr[l'^PL]} 

= tr[P EiLL"^}] 

= tr[P(0Q P"^ + AQ^^)] 

= n 0^ + tr(P A Q A^) ( 3 . 3 - 186 ) 

U Jv 

and 


E{X 2 U) = E{L^(AQ^'^ -1 

2 „-lv- 

- Oq P > 

^ AQ^A^PL} 

= tr[(AQ^^ + 

2 „-U-l 
Oq P ) 

AQj^A^P 

E{LL^}] 

= tr[(AQ^^ + 

2 ,-^-1 
“0 ^ > 

AQjjA^P 

(AQ^^ + Oq P"^] 

= tr[AQ^'^P] 




= tr[PAQj^^] 



( 3 . 3 - 187 ) 


we find that 

E{V 2 PV 2 + X 2 MX 2 > = " ( 3 . 3 - 188 ) 


This leads us to an unbiased estimate for the BLE variance of unit 
weight 


A^ A ^ 

.2 + ''2'®2 

°0 n 


( 3 . 3 - 189 ) 


which can be shown to hold also for Q positive semidef inite . As men- 

A. 

tinned above and as will be seen in the simulations of Chapter 4 , this 
a posteriori value will indicate the compatibility of the baseline 
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observations and the geophysical model introduced through the moment 

matrix Q . 

A 

The physical meaning of the BLE corresponds to one of the 
alternate approaches discussed in Section 2.4.3. In this case 
information is available on the expected deformations of the poly- 
hedron stations, for example, from an adopted absolute motion plate 
model (see Section 4.2). This would essentially eliminate the 
singularity of the CTS problem since the expected deformations refer 
to an absolute frame of reference fixed in the mantle (or crust and 
mantle) . In order to improve on these expected deformations and to 
test the deformation model, geodetic observations are taken (repeated 
baseline lengths) . By estimating the deformations of the polyhedron 
stations periodically, the reference frame is maintained since the 
CTS stations can now be assigned updated positions in that frame. 

Any adopted model should also include model parameter standard 
errors. This information and the expected deformations themselves 

are included in the moment matrix Q for the BLE. Another approach 

A 

is to simply correct the station coordinates directly and then use the 
stochastic portion of in estimating any residual motions that remain. 
The next estimator will follow this approach. 
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3.4 Unbiased Estimation 


3.4.1 Bayesian Estimation 

In the previous section, two possible biased estimators were 
presented. In the case of no a priori Information on the deformation 
vector X , these reduce to the ordinary pseudoinverse (or equiva- 
lently Inner constraint) estimate, X = n"*^U. In the case of the avail- 
ability of a Q matrix, we Investigate whether an unbiased estimate of 
X exists, and under what assumptions. 

Consider another extended GGM model (L, AX, Q^, X, Zjj) similar 
to the BLIMBE estimation model except that the moment matrix Q- of 
(3.3-37) Is split Into X and E-. For random variables X and V 

X. 

E{x} = X ; D[X] = Z- » E{(X-X)(X-X)'^} (3.4-1) 

X 

and 

e{v} = 0 ; D[V] = = E{w'^} = (3.4-2) 

Given the probability distributions of V and X we obtain the con- 
ditional distributions (Chlpman, 1964) 

E{l|x} = AX ; D[l|x] = E- (3.4-3) 

E{l|v} = ax + V ; D[lIv] = AZj^^ (3.4-4) 

It follows that the unconditional distribution of L Is 
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e{l} = E{e{l|v} = E{e{l|x}} = AX (3.4-5) 

D[L] = E{D[l|v]} + DfE{L|V>] 

„ (3.4-6) 

Compare these results to (3.3-141) and (3.3-142). 

In order to introduce the random expected deformation vector 
directly and not through the moment matrix as for BLIMBE and BLE we 
define a new random vector 

= X - X (3.4-7) 

There now can be written two sets of observation equations 


where 

and 


A linear estimator of X is given as a combination of L and L by 

A 

(3.4-11) 


X = GL + = [G G^] 


L 


A 


V 


= 


X + 




I 




(3.4-8) 


4 = ^ 


(3.4-9) 





V 


'Hy o' 

E^ 


► = 0 ; D 


= 





V,, 


0 Zz 


^XJ 


X 


X 


(3.4-10) 
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Let us construct an estimate that is unbiased and of minimum variance 


following (Chlpman, 1964). From (3.4-1) and (3.4-5) 

E(X) = GAX + G X 
A. 

= (GA + 

By equating to X , for unbiasedness 


(3.4-12) 


Gx = I - GA 


(3.4-13) 


The minimum variance condition then requires 


tr^[G I - GA] 




^X 


(I - GA) 


minimum (3.4-14a) 


or 


tr{GEj^G’^ + (I- GA)Z-(I-GA)’^} = minimum 


(3.4-14b) 


i.e., the trace of the covariance matrix is minimized. Comparing this 
expression to the risk function (3.3-176) for BLE it follows that 


G 




(3.4-15) 


This differs from G„ of (3.3-146) in that Z- replaces Q . It follows 

^ A A 

from (3.4-9, 11, 13, 15) that the estimate for X , call it X^, is 
given by 


X3 = z/(AZj^%r^)-^L 

+ [I-ZjjA’’(AZj^'‘^ + Zj^)~^A]X 


(3.4-16) 
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Thus, has the properties of unbiasedness, under the assumption 
E{X} = X, and minimum variance (but Is not BLUE being heterogeneous) . 

In other words. It has minimum mean square error In the class of unbiased 
estimators. Recall that the BLE has the property of minimum mean square 
error but In the class of biased (homogeneous) estimators. 

For positive definite E- 

M = (3.4-17) 

X^ = (N+M)"^U + [I - (N+M)"^N]X (3.4-18) 

= X* + [I - (N+M)~^N]X (3.4-19) 

/\* 

where Is similar to the BLE (3.3-150) but with E- Instead of Q^. 

/s* 

The second term In (3.4-19) Is the bias of X 2 . 

A 

Estimate X^ (In 3.4-18) can be written In another form as 

X^ = X + (N+M)“^A^P(L- AX) (3.4-20) 

which Involves a correction term to the expected deformation. In 
other words, the second term In (3.4-20) can be considered an estimate 
of the residual deformation 

X^ - (Xg + X) (3.4-21) 

where X^ are the CTS coordinates at an epoch t . Examination of 
(3.4-20) Indicates that linearization of the mathematical model 
(3.2-1) Is about the fundamental coordinates Xq. In the no noise case 

L = AX (3.4-22) 
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and 


= X (3.4-23) 

Finally, in yet another form, for positive definite and using 
identity C3. 3-149) on C3.4-16) 

= (N+M)"^U + [I- Ej^)"^A]X 

= (N+M)"^U + [m“^-m"^a'^(AM“^a'^+ Z,)~^AM”^]MX 

^ (3.4-24) 

= (N+M)"^U + (N+M)~^MX 
= (N+M)“^(U + M() 

the familiar weighted parameter estimate (Uotila, 1973) or the Bayesian 
estimate (Bossier, 1972). 

The covariance matrix is given by (3.4-14) and (3.4-15) as 

E- = GE + (I - GA)E-(I - GA)"^ 

A 3 L A 

= GEj^G*^ + ~ ®^X " 

= G(AE-a'^+E_)G^ + 2- - GAE- - E-aV 

A 1j a a a 

= E^^(AEjjA'^ + Ej^)"^AE- + E- (3.4-25) 

= ^x - V^'^x 

When E— is positive definite, using identity (3.3-149) , 
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(3.4-26) 


^X3 = 

Since the Bayesian estimate results from a combination of two 
tyiJes of measurements (see (3.4-8), associated with each type could be 
a different variance of unit weight. Define 

= E{VV*^} = QqR (=cTqP R positive definite) (3.4-27) 

E- = S positive definite) (3.4-28) 

(3.4-29) 


Then, estimate (3.4-16) may be written more generally as 

2 T ? T 2 —1 

X 3 = TQSA^(TQASA^ + aQR) L 

2 T a2 T 2 -1 — 

+ [I-TqSA^(TqASA^ + 0qR) ^A]X 

T T 2 -1 

= SA (ASaSp R) L 

+ [I - sa’’'^(asa’’+p^r)"^a]x V 
In the case of positive definite R and S (Chipman, 1964) 


(3.4-30) 


(3.4-31) 


/\ m 2 -1 T 

X 3 = [A PA + p M) A PL 

+ [I - (a'^pa+p^m)"^a'^pa]x 


(3.4-32) 


Different assumptions could be made about 0^ and Tq as outlined in 
(Bossier, 1972). For the purpose of this investigation we assume that 

Tq = % (3.4-33) 

2 2 

(the case when will be investigated in a future study) so that 

for the a posteriori variance of unit weight (Bossier, 1972) 
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(3.4-34) 



where , in general , 

/sT ^ ✓v'T /V 

[V PVl = V PV + 

= (L-AX)'^P(L- AX) + (L^-X)'‘^M(Ljj-X) 


(3.4-35) 


and 


/sT *r X /N'T at 

[V PV] = lJpl^ + - X U - X^ML^ (3.4-36) 

when M Is positive definite. In the expression for n is the 

number of geodetic (baseline length) observations. Note that for posi- 
tive semidefinite Z-, we use M = for (3.4-35). 

A 

The physical interpretation of the Bayesian estimate (X^) is 
the same as that of the BLE (X 2 ) explained in the previous section. 
They differ in how the expected deformations X are incorporated into 

A 

the adjustment. For X^, the expected deformations are added directly 
as corrections to the CTS coordinates (see (3.4-20)). Therefore, the 
estimated parameters are residual deformations. Since the deformation 

A 

model is so directly introduced we refer to X^ then as a "strong" 

A A 

Bayesian estimate. For X 2 (and X^ - BLIMBE) , the expected deformations 
are introduced in a weaker manner via the moment matrix Q and the 

A 

estimate is the total deformation. Therefore, we refer to X 2 as a 
"weak" Bayesian estimate. In fact, BLE is sometimes referred to as the 

A 

Bayes (instead of Best) Linear Estimate (Rao, 1976). Furthermore, X^ 
is an unbiased estimate though under the "strong" assumption E(X) = X, 
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while for the biased estimate X 2 we make no such assumption. Finally, 

/N ^ As. 

for X^, X is assumed deterministic and X stochastic as for X^. For 

As 

X 2 the opposite assumption is made. These differences will be studied 
in the simulations of Chapter 4. 


3.4.2 Beat Linear Conditionally Unbiased Estimate (BLICUE) 

As seen above, we have found an unbiased estimate for X under 
the assumption that E(X)= X. In this section, we present a condi- 
tionally unbiased estimate following essentially (Plackett, 1950; 
Chipman, 1964; Theil, 1971; Bossier, 1972). 

The starting point is the model (L,AX|CX = CX,Zjj»Qy) for the 
two sets of observation equations 


L 


A 


V 


= 


X + 


_^X_ 


C 




where C is the constraint matrix (3.3-90). The second set of equa- 
tions contain the weighted constraints 


CX = CX + Vjj ; (Ljj = CX) 


(3.4-38) 


so that 


Vj^=C(X-X) (3.4-39) 

where X is the expected deformation of the polyhedron stations. For 
this model 
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IV V 

^ = 0 ; D 


■ 2 -1 

0 cz-c^ 


(3.4-40) 


where 


D(V^) = E{V^vJJ} = CZ-c'^ 


(3.4-41) 


so that 


X ; D 




0 CZ-C- 


In contrast to the Bayesian estimate for which each station is 
treated individually in applying the expected deformations X, here 
the a priori information is reduced to six constraints. For example, 
the sum of changes in the X coordinates of the CTS stations may not 
sum to zero. Of course, if C = I then this model is equivalent to 
that of X^. 

The least squares minimum variance solution of (3.4-37) is 


X^ = [N + C^P^C]"^[U + C^Pj^CX] 


where 




Since AC = 0, it can be shown (Chipman, 1964) that using the notation 
of section 3.3.1 


+ T T 

= [N + C^P^C] Vp 
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and 


Cp I = [N + (3.4-46) 

X 

Recall that this means that (3.4-45) and (3.4-46) satisfy the four con- 
ditions of a minimum M-norm P-least squares g-inverse where in this 

case M = I. For (3.4-46), _ is P-least squares. Therefore, X 

Pxi X 

can be expressed as 


= A^jL + Cp pCX (3.4-47) 

X 

or in a more revealing manner as 

X. = n'^U + Cp CX (3.4-48) 

4 P^i 

y\ 

applying a result similar to (3.3-33). Thus, we see that X^ is com- 
posed of the BLIMBE with the condition CX « 0 (M = I) and a correction 
term that introduces the deviations of the reference frame conditions 
CX from zero due to the possible secular deformations of the poly- 
hedron. If no deformations are expected (X = 0) at an uncertainty 
level given by E-, X^ reduced to the standard pseudoinverse solution. 
Now, 


E{X^> = [N + c’^^P^C]"^ [NX + C^P^CX] 
= X 


(3.4-49) 


so that X, is an unbiased estimate conditional on 
4 


E(CX) = CX 


(3.4-50) 


Since it also has the minimum variance property (in this class of 
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estimators) we shall refer to as the Best Linear Conditionally 
Unbiased Estimate or BLICUE for short. Theil, (1971) refers to this as 
the mixed estimator and shows that it has a Bayesian inference inter- 

A 

pretation which is apparent when comparing it to X^. Therefore, it can 
be interpreted as a combination of the BLIMBE and Bayesian approaches. 
Furthermore, it bridges the gap between weighted parameter and con- 
strained estimation being a weighted constraint estimate. 

For the computation of the BLICUE by (3.4-43) only Cayley 

algebra is required as long as R(?y) = R(C) = 6 in which case 

T ^ 

N + C Pj^C has full rank. Otherwise, X^ could be computed from 

(3.4-48). 

The covariance matrix for BLICUE is derived as follows. Con- 
sider 


where 


= [N + c^(ci:-c'^)“^c]"^a’^p 
G 2 = [N + c'^(c2:-c'^)“^c]"V(cz-c^)"^ 


Then, 



where 


Z 


L 



(3.4-51) 

(3.4-52) 

(3.4-53) 

(3.4-54) 

(3.4-55) 
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(3.4-56) 


h = (cz^c ) 

Lx X 


which reduces to 


+ c'^(CPxC^)"^C] ^ 


(3.4-56) 


Deriving an expression for V PV yields an Interesting result . 


Denoting 


/v'T ^ A'T A A'T* A 

V PV = V^PV^ + VxPxVx 


(3.4-57) 


the second term is seen to be 

V^Px^x = [C(X-X)]'^(CZ-c'^)“^[C(X-X)] = 0 (3.4-58) 

since the weighted constraints require that CX = CX. Then, 

V PV = V,PV, = L PL - X, U (3.4-59) 

11 4 

which is the same as for the BLIMBE, and therefore, so is (3.3-56). 

By examining (3.4-48) one can interpret the BLICUE approach to 
monitoring deformations and maintaining the CTS, as dealing with the 
"systematic" part of the deformations in a P -least squares sense and 
the remaining "random" part in a P-least squares sense. As stated in 

A 

Chapter 2, applying the constraints of the form CX = L (L ^ 0) does 
not imply that the deformations Include global motions, rather secu- 
lar motions that do not average out. That is, for the expected defor- 
mations X, with respect to the reference frame implicit in the 
geophysical model, and for a particular station distribution, one 

expects that CX = L and not CX = 0. 

X 
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Table 1. Properties of Deformation Estimators 
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3.5 Summarizing the Properties of the Four Estimators 


Table 1 summarizes the properties of the four estimators pre- 
sented above. In the next chapter, numerical comparisons will be made. 

3.6 Addition and Temporary Deletion of CTS Stations 

3.6.1 Introduction 

The CTS frame is defined at an initial epoch by the adopted set 
of coordinates of a polyhedron of stations. In order to maintain the 
system, the CTS coordinates are updated periodically for the deforma- 
tions of the polyhedron. It is very possible that from time to time, 
one or more of the stations will not be able to participate in a parti- 
cular deformation analysis observing session. Furthermore, it must be 
anticipated that new stations will be added to the system periodically. 
Both of these occurrences must be accounted for, in order to maintain 
continuity and avoid ambiguity in the reference frame definition. We 
will adopt a least squares collocation approach to handle these situa- 
tions. 


3.6.2 Addition of New CTS stations 

In order to deal with this situation, we will apply the general 
model for least squares collocation (Moritz, 1980b) 

L = A^X + t + V (3.6-1) 

The observation vector L is composed of different parts. The first 
contains X , a non-random parameter vector to be estimated and A^ the 
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usual design matrix. The second is the signal part t, of a random nature, 
and third, a noise vector V due to errors in the observations (possibly 
also model errors). In our application, X yields estimates for the new 
CTS station coordinates. The signal portion is 

t = A^S (3.6-2) 


where the signals S are the deformations to be filtered, and their 
design matrix. The L vector is defined as before. Then, the expanded 
linearized equations for the mathematical model (3.2-1) are 


L = A^X + A 2 S + V 


= [A^ A2] 


+ V 


(3.6-3) 


In order to isolate the stochastic portion of (3.6-3), define 


L = L - A^X 


(3.6-4) 


The covariance function (moment matrix) for L is then 


Q- = E {Ll'^} 

= E {(A 2 S + V)ik^S + V)^} 

= A 2 E {SS^} A 2 + E{W^} (3.6-5) 


where we assume that signal and noise are uncorrelated. It follows 
that 

Ql = A2QgA2 + (3.6-6) 

where 

Qg = E {SS^} ; Qy = E {VV^} (3.6-7) 


are the signal and noise covariance functions (moment matrices) 
respectively. 


Ill 


Note that Q is constructed only for those stations whose deformations 
are being estimated while refers to all observations. The cross 
covariance function for the deformations and observations L is 




= e{s(a,s + V)^} 


= E{Ss’^}A2 

T 

= V2 


(3.6-8) 


The collocation estimates for X and S are then given by 
(Moritz, 1980b) 


X = [A][(A2QgA^ + 


(3.6-9) 

(3.6-10) 


and 


^X ^^1^^2^S^2 


(3.6-11) 




(3.6-12) 




(3.6-13) 


The design matrix 
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A 

n u 



n q n u - q 


has the general pattern 



Fig. 2 A-Matrlx Structure for Addition of CIS Stations 

The number of rows in I is the number of baselines in the original 
polyhedron. The number of rows in II is the number of baselines with an 
original station at one end and a new one at the other. The number of 
rows in III is the number of altogether new baselines . Since both A^ 
and A^ are computed from the same baseline length model, they are both 
rank deficient. This explains the pseudoinverse in the expression for 
X (3.6-9). 
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Reiterating, the vector X contains the estimated coordinates 
of the new CTS stations. The linearization, in this case, is taken 
about some approximate station coordinates. The vector S contains 
the filtered signals (deformations) relative to the expected deformation 
and the linearization is taken about X^, the fundamental coordinates. 

It should be mentioned that once new stations have been added to the 
system, in subsequent deformation analyses the linearization of the 
mathematical model (3.2-1) should be taken about X^ where t is the 
epoch at which the new stations were added. In any case, it may be 
useful to update the linearization point periodically from convergence 
considerations although the deformations are small compared to the 
baseline lengths. 

An examination of the estimates for X and S is quite 
revealing. Consider for a moment that there are no new CTS stations. 
Then 

S = QgA’^(AQgA'^ + Q^)"^L (3.6-14) 

which is just the BLE estimate (3.3-146). Thus, application of this 
estimator represents the philosophy of considering the deformations as 
signals to be filtered from the observational noise. Recall that in 
the estimation model for BLE, the deformation vector X was assumed to 
be random. Next, consider the X vector when there are no signals 

X = (A^Q'^A)'^ A^ q"^L (3.6-15) 

Adding an M-norm in the parameter space yields the BLIMBE (3.3-43 ). 
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Thus, the BLIMBE approach follows the philosophy of considering the 
deformations as parameters, i.e. deterministic quantities. 

3.6.3 Temporary Deletion of Several CTS Stations 

Some CTS stations may be unable to participate in a particular 
deformation analysis observing campaign. Provisions must be made for 
this in order to maintain the reference system. In this case, it is 
possible to apply the prediction capabilities of collocation, which is 
equivalent as shoim by (Dermanis, 1976) to minimum mean square error 
prediction. Recall that the BLE yields minimum mean square 
error so that this method will be an extension of the BLE model. 

Consider the BLE model with new notation reflecting its fil- 
tering interpretation 

L = Bt + V (3.6-17) 

where now, t represents the signals (deformations) actually measured 
(at the participating CTS sites) and B their corresponding design 
matrix. Then 

C. = B Q^b"^ + Q„ (3.6-18) 

u t V 

from which 

S = QgB’^(B Qj.b'^ + Q^)”^L (3.6-20) 

Here S , the signal vector has two components 
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(3.6-21) 


S 



where t are the deformations to be filtered at the participating 
stations and u are the deformations to be predicted at the missing 
sites. The structure of the B matrix Is as follows 



>- m q 
m 


0 

m u-q 


where there are m baseline observations. The signal covariance matrix 
Is constructed from the geophysical deformation model for the q 
observing stations. The full Q matrix Is computed for all the CTS 
stations whether they have observed or not. Actually, both Q and Q 
are subsets of the global Q matrix that can be computed for any point 

on earth (compare to M of section 2.4.3). 

G 

The prediction Is accomplished through the adopted geophysical 
model from which Is derived the signal covariance function. In this 
context the covariance function has probabilistic justification. 

Dermanls (1976) shows that the choice of Inner product (l.e. weighted 
norm) Is equivalent to the choice of covariance function In colloca- 
tion. 


3.7 On the Estlmablllty of the Baseline Length Change Estimates 

The question raised In this section Is whether or not the 
"adjusted" baseline length change vector 

L = AX (3.7-1) 
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is estimable, or in other words is it an unbiased estimate of L for 
each of the four estimators? In our context, is the estimated change in 
the size and shape of the polyhedron unbiased? 

For BLIMBE, 


= AX^ = AQ-N(NQ-N)®a’’pL 
= AG^L = AaJ^L 


(3.7-2) 


Then, 


E{£^} = AA^^{L} 


(3.7-3) 


Therefore, is an unbiased estimate for L though is biased. 

For MINDLESS (differentiated from the BLIMBE by an asterisk) 

£* = AX* = A(N+Q|)“^N[N(N-fQ|)“^N]®A'^P(AX + V) 

= AN%(N+<5|)"^N[N(MfQ|N]®[NX + a'^PV] (3.7-4) 

= AN%X + AG*V = AX + AG*V 

using identities (3.3-114) and (3.3-115). It follows that 

E{£*> = ax = L (3.7-5) 

Therefore, Lj^ is unbiased, too. 

Recall that the MINOLESS and BLIMBE are equivalent when Q- is 

* 

positive definite. Note that AG^ and AG^ are idempotent and, therefore. 
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both are projection operators. That is, they map the observation vec- 

/V 

tor L into L an element of the column space of A . 

For BLE, 


/s ^ T T 2 -1 -1 

L2 = AX2 = AQ^'-CAQ^^ + 

= AG2L 


(3.7-6) 


and 


E{AX 2 > = AG 2 AX ^ AX. (3.7-7) 

^ A 

in general. Therefore L 2 is biased as is X 2 . Note that since AG 2 is 
not idempotent, it is not a projection operator. 

For the Bayesian estimate, 

£3 = AX^ = AX + AZ^^(AZ^^ + aQP"^)”^(L- AX) (3.7-8) 

and 


E{L^} = E{AX^} = ax = L (3.7-9) 

if E{X} = X, the assumption made for this estimate. Similarly, the 
BLICUE 


E{L^> = L (3.7-10) 

under the assumption E{CX} = CX. Therefore, both estimates yield con- 
ditionally unbiased estimates for X and L . Note that both esti- 

/V /V 

mates are not homogeneous and, thus, one cannot express X and L as 
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X = GL or L = AGL, as can be done for BLIMBE and BLE, except by adding 
another term to L (see (3.4-11) and (3.4-43)). 

How does one interpret the fact that is unbiased and 
biased? For BLIMBE, it is apparent that the model moment matrix Q- has 
no effect on the estimation of L . That is, is Invariant with 
respect to a weighted norm in the parameter space and would be the same 
as obtained in a free adjustment without any geophysical information at 
all. This also follows from examining the BLIMBE estimation model 
(3.3-35)- (3.3-40) in which there is no connection between L and X . 
On the other hand, for BLE, L 2 is influenced by the a priori information 
given by Q , or in other words by the expected value L of the new 
distances computed from X the expected deformations. Therefore, it 
comes down again to what is preferred a biased or unbiased estimate for 
some parameter, or does one prefer to ignore a priori information or 
not. 

It is well known that a biased estimator can improve upon 
unbiased estimators if there is some a priori information about the 

/s. 

unknown parameters (Rao, 1973). In our application we can show that L 2 
has minimum mean square error. That is, 

MSE(L 2 ) = E{(L 2 - L)(L 2 - L)'^} = minimum (3.7-11) 

From (3.7-6) 
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(3.7-12) 


MSE(L2) = E{[AG2L-L][AG2L-L]'^} 

= E{[AG2~ I]Ll'^[AG2" I]^} 

= E{[AG2-I][AX + V][AX + V]'^[A G2-I]'^} 

(using Qyy = 0) 

= [AG2 - I]AE{Xx’^}a'^[AG2- I]"^ + AG2E{w'’^}G2a’’‘^ 

= A[G2A- I]E{XX^}[G2A- u'^a'^ + AG2E{W^}G2a’’^ 

= AE{(X-X)(X-X)'^}a'^ (3.7-13) 

= A MSE(X) a"^ (3.7-14) 

which follows from (3.3-128). Since the BLE estimator minimizes the 
mean square error (in the class of biased estimators ) , then it follows 

A 

directly that (3.7-11) holds and improves upon L^. Using the same 

A 

reasoning it follows that for the strong Bayesian estimate has mini- 
mum mean square error (also minimum variance - being conditionally 

unbiased) in the class of conditionally unbiased (heterogeneous) 
estimators.- 
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4 DEFORMATION ANALYSIS SIMULATIONS 


4.1 Introduction 

In Chapter 3, four possible estimators were presented as pos- 
sible candidates for CTS deformation analysis. Each one has its parti- 
cular optimal properties. In this chapter, we study the significance 
of these properties by running a series of simulations. The main ques- 
tion to be answered is whether or not an absolute motion plate model 
should be adopted, as opposed to performing an ordinary (unweighted) 
pseudoinverse solution (a free adjustment). In addition, we seek to 
determine the best way to combine geodetic baseline observations and 
geophysical models in the estimation of crustal deformations. Recall 
that in a free adjustment, the singular normal matrix is augmented with 
the constraint matrix C such that CX = 0 (alternatively, pseudoinverse 
algebra is used to solve the set of normal equations. That is, these 
constraints arbitrarily impose no net translation or rotation for the 
estimated deformations X without any real physical justification. 

/V 

The introduction of a geophysical model can direct X to a physically 
more meaningful solution. A similar concept has been applied in 
(Prescott, 1981) for monitoring deformations along a strike slip fault. 

We "adopt" the absolute motion model AMl-2 of (Minster and 
Jordan, 1978). Since using this model leads to a positive semidefinite 
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moment matrix Q , we compare the four estimators first using a positive 
A 

definite and then a positive semldeflnlte model matrix. 

4.2 The Model Matrix 

As mentioned in section 2.4.3, it is postulated by Wilson and 
Morgan that hot spots (or sources of volcanic magma) are fixed in the 
deep mantle and, thus, form a set of rigid points that serve as an 
absolute reference frame for plate motion. In AMl-2, rotation vectors 
for eleven major plates are given in this absolute frame. The compo- 
nents of these vectors are the poles of rotation of the plates ((J>p,Xp) 
and their rates of rotation (o)p) . Standard errors are also provided 
for these estimates in Table 2, for which the data is obtained from 
(Minster and Jordan, 1978). As can be seen, the errors in the 
slow moving plates are quite large for the pole parameters. The 
intention is not to advocate a particular model but to adopt one for 
simulation purposes. We note that there is controversy regarding the 
hot spot hypothesis (e.g., Le Pichon, et al. ,1973). Nevertheless, 
as mentioned before, Bender (1981) points out that absolute motion 
models derived from different geophysical assumptions about plate 
motions differ by about 1 cm/year. 

A model moment matrix for the expected absolute velocities of 
the CTS stations can be constructed from an absolute motion model as 
follows. The velocity of station i on plate j in an absolute frame 
is given by (Minster and Jordan, 1974) 

X (4.2-1) 
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MillioQ Years 

Not used in the simulations 



where 




COS(|)jCOSXj 

cos(|)j8lnXj 


slruj). 


(4.2-2) 


is the j'th plate's rotation vector, and 




R 


cos <|)^COS X^ 
cos<j)^sinX^ 
sin((). 


(4.2-3) 


is the spatial Cartesian I'th station vector. A spherical earth 
approximation is sufficient for the description of plate kinematics 
where R is the radius. Then, 


V.. = 
ij 


V 


"ij 


V 

= Ru). 

^Ij 

3 

V 


L "ijj 



cos<j)jSin4i^slnXj - sin(J>jCos())^sinX^ 
sin((>jCos4)^cosX^ - cos(j)^sln4)^cosXj 
cos<{)jCos({)^sin(Xj^ ~ 


(4.2-4) 


By error propagation, the variance-covariance matrix for the 
station velocities is given by 


where G is the partial derivative matrix of the form 


(4.2-5) 
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(4.2-6) 







M 


) 


i = 1, 

j = l, 


,k 




where k denotes the number of stations, Z the number of plates and 
Zp Is the covariance matrix for the plate parameters given from 
Table 4.1. Here we assume a diagonal covariance matrix. It has been 
learned that a full covariance matrix is available (J. B. Minster, 1982, 
private communication) although too late to be used for these simula- 
tions. The partial derivatives are computed as 


8V 


3<|). 


Eii). [-sin<j).sin(t).sinX, - cosd). cosA.sinX. ] 
J J 1 J J i 1 


3V 

Jhi 

9X. 


Buij cos({)j sin((»^cosX^ 


8V 

X. . 

11 

9w. 

J 


R[cos<|)jSin4)^sinXj - sin(J)^cos(j)^sinX^] 


3V 




l±i = 


E(jj. [cosd). cosd), cosX, + sincj) . sind» . cosX , ] 
j ^i i ^3 ^i j 


3V, 


— ^ = Ro),cosd).sind).sinX. 

J 3 


3V, 


3 ( 0 , 


-±1 = 


R[sind>jCos({)^cosX^ - cos(J)jSin(j)^cosXj ] 


9V, 


3d>, 


11= - 


R( 0 j sin(})j cos(()^sin (X^, “ ) 


(4.2-7) 

(4.2-8) 

(4.2-9) 

(4.2-10) 

(4.2-11) 

(4.2-12) 

(4.2-13) 
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(4.2-14) 


9Vz 

= -Ru) cos(J) . cos<{) cos (X . - X , ) 
oA. j 1 i 1 j 


9V. 


Jl = 


Rcos(() cos(|) sin (X. -X.) 

oWj 3 1 1 j 


(4.2-15) 


It is apparent that 



if station 1 is not located on plate j . 

We assume the deformation vector X to be related linearly to 
the velocity vector by 


X=V(t-tQ) (4.2-16) 

where t - tQ is the time elapsed from the initial CTS epoch. The model 
moment matrix for the deformations is given by 


= E{XX'^} = (t- tQ)^E{w'^} 
= 

using 

= E{w’’‘^} = 


(4.2-17) 


(4.2-18) 


where V denotes the expected velocity vector computed from (4.2-4), 

X is the expected deformation vector computed from (4.2-16). and and 
are the covariance matrices of V and X , respectively. Recall 
that in Chapter 3 the moment matrix Q is introduced differently for 

A. 

each estimator. Note that in (4.2-16) and (4.2-17) as the time inter- 
val gets shorter the deformation gets smaller and so does its 
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covariances. This is somewhat misleading since the plate models are 
given as long-term average rates (over approximately 50 years to several 
million years (J. D. Minster , 1982 private communication)). However, 
Bender (1981) suggests that the present rates of motion should not be much 
different from the long-term average rates. This is then the assumption 
made in these simulations and the one that will eventually be tested by 
geodetic observations. 

The matrix for a horizontal plate model, such as AliL-2, is 
always rank deficient. Implicit in the plate velocities (4.2-4) are 
horizontal motions between rigid plates. That is, no vertical motion is 
indicated and baseline lengths on the same plate should not change. The 
matrix is given though in terms of deformations in 3 components 
(x,y,z) per station. If we would express the deformations in a local 
system (see section 2.4.1), the height component would drop out. There- 
fore, there is one rank deficiency per station due to the vertical com- 
ponent. In addition, for stations on one rigid plate there is one 
rank deficiency per non-redundant baseline. For example, consider 4 
stations on one plate. The model predicts 12 deformations, 3 per sta- 
tion. For 4 stations, a quadrilateral, there are five independent 
baselines out of six. Therefore, there are 9 rank deficiencies — 

4 vertical motions + 5 rigid independent baselines. This reasoning 

X 

applies to the Z portion of Q . The XX portion always has a rank of 
X X 

one so that the sum of the two matrices, i.e., Q , can at most be 

X 

increased by one over the rank of Z . 

X 
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The above example of 4 stations Indicates also that the minimum 
number of stations per plate should not be less than four from the 
reliability point of view. For n < 4 there is no redundant baselines 
and therefore no check for systematic errors or site stability prob- 
lems . 


4.3 Numerical Tests 

4.3.1 Simulation Procedure 

In Chapter 3, four estimators were presented as possible candi- 
dates for deformation analysis. Each one has its optimal properties. 

In order to get a better feeling for the applicability of these esti- 
mates, this section describes some numerical tests. For example, it 
was shown that the BLIMBE is a minimum bias estimator. However, if this 
minimum bias is large, it makes little sense to use this estimate for 
deformation analysis. Through controlled experiments described below, 
it is possible to assess the magnitude of the bias, and similarly test 
the other properties of the different estimators. 

Several numerical comparisons were made for the four estima- 
tors. The least squares property is examined through the computation 
of V PV and the minimum norm property by X MX. For the two biased esti- 
mators, the trace of the bias matrix is computed, l.e., 

tr[(I-GA)Q^(I-GA)'*^] (4.3-1) 

where for BLIMBE (Q positive definite) 

A 
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(4.3-2) 


(GA)i = QjjN(NQ^N)‘^N 

For the minimum M-seminorm P-least squares estimate (MINDLESS) with 
positive semidefinite Q 

(GA)^ = (N+Q+)"\t[n(n+Q+)"^N]‘^N (4.3-3) 

For the BLE, Q positive definite 

(GA)^ = (N-H)"^)"^N (4.3-4) 

and 

(GA)2 = Qj^'^(AQj^'^ + Z^)"^A (4.3-5) 

for positive semidefinite Q^. The minimum variance (mean square error) 
property is reflected by the trace of the different covariance (or mean 
square error) matrices. For the BLIMBE this quantity is added to the trace 
of the bias matrix to yield the mean square error. The ratio of maximum 
and minimum eigenvalues (referred to as the C-measure in section 4.3.4) 
assesses the condition of the covariance or mean square error matrices . 

Another property, and a most essential one, is how close does 
the estimate come to the true value. This is tested in a simulation 
environment where the "true" value is known. For this purpose, we com- 
pute the root mean square of the deviations of the estimated value from 
the true value as 
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(4.3-6) 


(4.3-7) 


where the bar denotes an average and p is the number of stations. Both 
variance and bias terms are computed and added to yield what we will 
call the RSMPE (the root mean square of the sampling error) . 

The simulation procedure is as follows. For a particular net- 
work at an initial epoch t^, we assign a set of "fundamental" coordi- 
nates Xq and compute their corresponding baseline lengths. Next, by 
means of the AMl-2 absolute plate motion model, we compute the expected 
deformations X = (X^ - X^) at some later epoch and compute the corres- 
ponding expected baseline lengths. To these v/e add 3 cm Gaussian noise 
and subtract from them the initial baseline lengths to yield the 
"observed" baseline length changes. We assume that the baseline 
lengths are re-observed after two years (t - t^ = 2 years). 

The adjustment algorithms outlined in Chapter 3 are followed. 

In all cases, the linearization of the baseline length mathematical 
model (3.2.1) is taken about Xq. It is assumed that the re-observed 
baseline lengths are accurate to 3 cm and uncorrelated. This is a 
reasonable assumption considering that at the present, it is possible 
to estimate individual intercontinental VLBI baselines with nearly such 
precision (Herring, 1981). The model moment matrix is computed on 
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the basis of the AMl-2 model as explained in the previous section. 

The network depicted in Fig. 3 is the basis for the simula- 
tions. Three stations were chosen per each of six major tectonic 
plates (North American, South American, African, Eurasian, 
Indian-Australian, and Pacific) and one station for two smaller plates 
(Arabian and Nazca) . The sites were chosen on the basis of several 
criteria. First, the station if possible should be operational or at 
least have been mentioned as a likely candidate. This criterion is met 
by at least eleven of the sites (stations 4-11, 13, 16, 17 of 
Table 3) . The remaining sites were chosen in stable parts of each 
plate according to Lowman's global plate tectonic map reproduced in 
Fig. 4, Furthermore, each plate should be well represented which was 
determined by an examination of the expected plate motion vectors shown 
in Fig. 3 and listed in Table 3. This same network will be used in 
the experiments of section 4.3.4. 

In order to compare the four estimators using first a positive 
definite Q matrix, a subset of 8 stations is chosen, one per each of 
the above mentioned tectonic plates (Table 3) . In order to construct 
a non-singular model matrix a secular vertical motion model was 
selected arbitrarily (and therefore not reproduced here) and its cor- 
responding model matrix added to that propagated from the AMl-2 model. 
It consists of vertical deformations in the range of ±3 cm per year with 
1cm accuracy. This eliminates the rank deficiencies due to the undefined 
vertical components. Since there is only one station per plate, there 
are no rigid baselines. The 8-station experiments are solely performed 


131 





'O 

(U 

4J 

0 

u 

(N 

1 


T3 

c 

cd 

M 

u 

o 

s 

4-1 

0) 

z 

c 

o 

•H 

4J 

CO 

rH 

3 CO 
fi 

•H 0) 
CO H 


0) 

cO 


H 

CO 


iH 

(U 

PM 

Q> 

1 

CO 

00 

'w' 

c 

CO 

0 

0) 

•H 

•H 

U 

•u 

CO 

•H 

4J 

O 

CO 

O 

1 

rH 

o 

0) 

CM 

> 

CO 


• 


00 


•H 





A 


cc 

> 

g: 

u 

o 


132 




133 


w 

0) 

•H 

4J 

•H 

U 

o 

r-( 

0) 

> 


tN 

I 


& 

c« 

T3 

a 

Jai 

U 

o 

a 

4J 

0) 

ai 

a 

o 

•H 


a 

Q 


•H 

Vi 

<V 

4J 

rtJ 

Ot 

I 

00 

(3 

0 

•H 

V 

■P 

cn 

1 

o 


n 

(P 

A 

«5 

H 



a 

•rl 


d) 

V 

rtJ 

04 


a- 
•a 
a 

4-)-pe 

oj 

O'rt O' 
OWO) 
o o 

►J 


(U C3 
t3x;H 
0^23 
*> U 
O O' 
•4JSB <U 

(0 a 

>-i 


a 

o 

•H 

+> 

tn 


t 

o 

as 


o »“ <N ^ ^ O' ^ cj- O ' 00 .o oo I vO ■■'n r». o o 

oo'0»-o<NLn»-LnLno-r>oo<Nirsj 

»— T— oooooT-oofor^oooiri'^^'^'N 

III I I I I 

O' (N 4” O irt f~- O' lO O' ^ or O >0 O O' tr» so 'JO 
O' on .'n ro <»' rsl »- O P'1 in f'J O' (N O' (N iD ro 

O 0000000*“ CO P*0'^^»“»~»~'~'“0 

I I I I I I I III III 

O' O' O' ('n O' ® pn r*- o so o- m po o If) P'1 

O iT O *” r" O Pn a' Psl O *“ P'1 •IT P*0 P'1 O' o ^ 

O O O O O O ^ P4 P'1 »” vO pj P'1 P'1 P'1 PP) ^ O IJS *“ 
I I I I I I I I I I t I I I I I I 


i uuu<*<<!*:as3uuu*:s:*;HHHO» 
04 Q) 04 « « 03 •< «« «a Du Hu (t4 >4 •« «< O O Q N < 
I Hufc4friOOO0OOUUUOOOS5SBas«J!05 

i< <« E4 H CO as as SB 04 cu 04 cn to c/) )-i H H as < 


p\j in 03 o oif) o pn o in o o *” so o ® «“ •“ O' o 
o*“m sr»" posrpp) r^poar pvimpnar 

CO «- po ro ® •- P'1 o CO pp) o po r- pp) O' m fN o vo 
psipn *~mf'i*-in®or»'0»“00'ir*"r-ino’ 
PP)^fNlPlP'lP'lP'U-PnPnP'|r-'r. (S 


o pp) o o in o o so o o o pn ^ in ® pp) vo in O' 

»— (M cs’inpp'pppp'pp) ppipppp)*” m tn 

O' o sO p* po o 3" o P '1 o p* pn o in O' ® p» ^ 


ps|pp| 

mm®®®=rp'i 

*- P'iPP)T-pnp'|. 

i-PsKN 

0 ) a 4J 

1 


1 II II 

1 

CO (TS 
0#»^-4H 
O 4 dcx* 

* 


WJ 


0 1 

# 


V 

« 

0 H® 

O' 


rH 0 ) 

# 

U 

M 

a « 

03# 0) 

<V 

<]; ITS d 

4) 

(TS ♦ 

H# Lj 

4-4 

XT3 a 0 

J3 

® w w 

O-r^ # 0) 

</l 

#> 9-H.H 

(fl 

♦ •MJj'r^'TS 

r^p-^«; « (D 

H# 

•HU x+j 

0) 

# H (TS d > u 

H 9 tn>-lt3 

« 

O H 0 flS 

d 

itSp^JS (TS (TS 0# 

•r-l ITS (4 W d (TS ITS 

>*M d 

0 OTU4J 


a O W<-l 0) OU3Q>l-l* 4JJ304 o O U CT'iO Q>-H 
<«MO(4l4SU 4 J-'^’p 1U) a rO O >4X)4J-^ 
43 -kI 0'(/]t 3 n)-H I (/} 9-J3 UlOa>UMUaU><4 
o (0 (4 o o>0 <04J ai(0i4!09aiohn)O(4u 
o u >-» o pj OT fc 4 fc< a* e H e tn ® u o >H ® w «4 


p" p< ® in ® r*» ® ® o *“ P'1 ® ar m ® r* ® ® o 

*■ *— r P'1 


JS 
o Jj 
O tJ 
O' (0 
Q 
■J1 

OrH 
U ITS 
PS U 

».r4 

o u 

4J 0) 

x: 

» 04 

w 

% 

0) • 
ITS 0) 
S5-H 
4J O 

■r4 a. 

O'. 

a a 
o o 

rH*H 
■iJ 
o iS 

o*J 
O 
U)U 
T3 
Ul U) 

«s» 

»JS 

o*> 

M 
ITS 
X Q) 

Ui 

w® 
*> u 
a (ts 


X 

u 

o 

a 

+1 

0) 

SB 

Q) 


r-l (3 Oi/i 
<0 0 0,1 
>4-4 ITS® 

* * 
« 


134 



for the sake of comparison with a positive semidefinite Q matrix. 

X. 

Next, the same set of experiments are performed for the 18-station net- 
work under realistic assumptions. Of course, the 18 station network is 
much stronger having more degrees of freedom. Using (3.3-60) the 8 
station network has 10 degrees of freedom, the 18 station network, 105. 

In order to study the effects of errors in the model matrix Q , 

X 

first 3 cm and then 6 cm Gaussian noise is added to the expected defor- 
mations (4.2-16) computed from AMl-2 in order to construct a weak 
but somewhat realistic model matrix. The simulated baseline length 
changes are computed as described above using the "correct" model. Thus, 
the geodetic observations detect the "true" deformations within 3 cm 
observational noise but the geophysical model is rendered somewhat 
incorrect, and inconsistent with the geodetic data. The noise level on 
the deformations was chosen according to the uncertainties attached to 
the AMl-2 plate rotation vector parameters. The propagated deforma- 
tions have standard errors on the order of several centimeters for 
t - tg = 2 years. This seems at first glance surprising considering 
the large uncertainties in the AMl-2 parameters. A closer examination 
indicates that the poles of the plate rotation vectors contain the 
largest uncertainties, particularly for the slower moving plates. 

These are not as critical as the plate rotation rates whose standard 
errors are comparatively smaller. 

4.3.2 Results of the 8-Station 8-Plate Experiments 

The results of the 8-station network experiments are listed in 
Table 4. it should be noted that all numbers are means over 3 runs 
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Table 4 Results of 8-Station 8-Plate Network Simulations (AMl-2 + vertical 
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* Model noise (cm) 

*• Nuabers in parentheses are root mean square values (cm) 
A Priori values 




each with different noise in order to get a more representative sample. 

More "noise loops" were not considered necessary since the results from 

the 3 runs did not appear to differ significantly. 

The first conclusion from these runs is that the results are in 

accordance with the optimal properties of each estimate as outlined in 

aT , 

Chapter 3 and summarized in Table 1. The V PV’s 'for BLIMBE and 

BLICUE are equivalent and minimum, both being a generalization of the 

ordinary least squares estimate (BLUE). It is interesting to note that 

T 

an error in the model matrix has a very small effect on V PV for both 3 
and 6 cm model noise. This error is reflected, as expected, in the norm 

A A<J> A 

X MX. This results from V PV being defined in the observation space 

A^ A 

while X MX is defined in the parameter space. As can be seen, the para- 

meter space variance of unit weight Tq (3.3-61 ) is a good indicator of 

the compatibility between the adopted model and the geodetic observa- 

''2 

tlons. The same holds for (3.3-189) for the BLE, where now both 

A*p A Alji A 

V PV and X MX are defined in the observation space. Furthermore, we 

A^ A A^ A 

recall that the BLE minimizes V PV + X MX, as is the case here. For 
the Bayesian estimate, the situation is similar to that of the BLE. 

A*JI A 

V PV is computed as in (3.4-36) and the variance of unit weight 
according to (3.4-34). Note that in Table 4 the contributions from 
both terms of (3.4-35) are listed. 

A^ A 

For the Bayesian case, X MX is not very informative. Recall 
that here M includes only the stochastic portion of AMl-2, as is the 
case for the BLICUE. In fact, for the latter none of the indicators in 
the table seem to reflect an error in the model matrix. 
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We see that this is approximately 0.9 cm for BLIMBE and 1.2 cm for the 
BLE. Note that It does not depend on the actual observations so that 
it is an a priori indicator and a function of the network geometry (the 
design), the observational accuracy and the geophysical model. In any 
case, though BLIMBE is minimum bias, that of the BLE is not much 
larger. The BLIMBE minimizes the bias at the expense of variance com- 
pared to the BLE. In other words, BLIMBE minimizes the variance in 
the class of minimum bias estimators while the BLE does so in the class 
of all biased homogeneous estimators. Furthermore, the BLE minimizes 
the mean square error in this class. However, the heterogeneous Bayesian 
estimator has even lower variance and thus mean square error (there is no 
bias) than both biased homogeneous estimators. 

The conditioning of the various solutions is reflected in the 
C-measure. In this example, the Bayesian estimator yields the most 
stable covariance matrix as can be seen from the correlation distribu- 
tion given in Table 5. As can be seen the BLICUE seems to possess a 
structure that lies between the Bayesian type estimators (BLE and 
Bayesian) and the pseudoinverse one (BLIMBE) as was indicated in sec- 
tion 3.4.2. 

Finally, we examine what may be the most interesting indicator, 
the RSl^E of the estimated deformation compared to the "true" 
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Table 5 Correlation Distribution for Deformation Parameters 
8-Station 8-Plate Networit 
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deformation. The magnitudes of the "true" deformations (computed from 
the AML-2 and the vertical motion model) for a two year period are 
listed in Table 6. Note the root mean square deformation of 6.3 cm. 

In all cases the bias term is small compared to the variance term. Of 
course, the Bayesian estimate yields the smallest RSMPE in the case of 
no model errors, i.e., when the assumption E(X) = X (3.4-1) holds. In 
this case, the adjustment Just involves filtering out the observation 
errors after the correct deformations have been applied directly 
(3.4-20). Note that the weak Bayesian approach of the BLE yields com- 
parable results even though the expected deformations are entered 
indirectly through the moment matrix Q . However, as the model errors 

iv 

increase, X moves away from X , the Bayesian estimate yields the 
largest RSMPE while the BLE is less affected. This will become more 
pronounced in the results of the next section. Here, the biased esti- 
mates and the BLICUE yield somewhat better results, i.e., they are less 
affected by model errors. 

4.3.3 Results of the 18-Station 6-Plate Experiments 

In this set of simulations, 18 stations are distributed, 3 per 
each of six major tectonic plates. The geophysical model is AML-2, 
thus the model matrix is positive semldefinite. The same series of 
tests were performed as for the 8-station network. This simulation is 
more representative from the point of view of greater redundancy (105 
degrees of function versus 10 for the 8-station simulation) and more 
interesting since only the AML-2 model has been used as conceivably would 
be done in practice. The results are listed in Table 7. 


140 


o 






O 

9VP 

n 




P 

0^ 




lOU 

OH 




>9 9 

•H 

a 

MOOfnovOMOOO 


MO 



*-MMini»OOOM 


9 


'M# 



9H 

■ « 

M 



09 


M'<* 1 

1 1 1 


PU 

oca 

O 


•H 

tMH 

P 



MM 

o 

U 



0 

QO 

« 

^nO•-<^Ovfr^MM 

a 

M3 

m 


flO»"(OOr*0*»M 

u 

0 CM 

TJH 

« 

tiftttti 


<9 W 

9*i 

AM 

MO*-m(ft»r>nM 

r* 

9 

+> 

o— 

'-I II 

M 

P » 

9 U 

•H 


1 

HO 

a9 

p 


VO 

9H 

9« 

9 



90 

o 

a 

^os^-ovroovo*” 

• • 

O Ou 

u 

M^ 

v009,90VMOO<^ 


pH 

<N 

0# 

tffftfl* 

o 

9 


IMM 

*»M9vOr»OMU> 

O 

o o 

0> 

O'- 

1 1 1 1 

a>H 

OH 

>9JM 

a 

1 

OP 

P 

ou 


9 

(Q9 

EO 



<na 

>9P 

» 



(PM 

MO 

H<4J 



*0 

9M 

9 O 



OP 


ow 

« 


o 

OU) 

H 

p 

U>9BUBHU0Q 

• tO 


+>® 

9 



J0 

M-P 

H 

PkSOUOBBPoa 

09 

MP 

O 9 

ca< 

PHaOt(OHBS<9 

OM 

-•M 

>H 



•H9 

9 

Ot 



P 9 

too 

>a 1 



99’ 

P 

ei(B 



a(0 

9(0 

9 



M 

O’® 

C3 


O* 

09 

9M 

MO 


M • 

P9 

09 

1 -H 

9 

O H 

(DO 

CM» 

»-p 

O 

P (0 O 0} 

COB 

ao 

S 9 

•M 

0 H H W 


OP-^ 

<«p 

P 

9 > 9H 

9P 

U 9 

tfj 

9 

A99 99M9 

90 

NO 

'Ml 

P 

9HO CM U «H 

90 

9 H 

OCD 

(0 

99 H OPP 

B« 

O %P 



jam (oouna 


HO 9 

SM 


O flp 9 9 M 9 M 


P '9 a 

90 


DOIk«B(OOM«e 


9 9H 

(O'M 

• 



aPM 

MHO 
0 9»M 

VO 

o 

^ GO oro vO ov o 


P90i 


9B 

^T— y r»fSj 


OO CM 

O 




OH 9 

pH 









« 

9 





H 






141 



Table 7 Results of 18-Statlon 6-Plate Network Simulations (AMI- 
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nodel noise (cm) 

A Priori values 

Nunbers in parentheses are root mean s-.juare values (cm) 



The first innnedlate conclusion is that the estimator of section 


3. 3. 1.4 for seminorms is no longer minimum bias. Recall that the proof 
of the minimum bias property assumedj a positive definite Q- matrix. 
Therefore, the first estimator in Table 7 is denoted as MINDLESS, in 

this case <a minimum seminorm P least squares estimate. Recall that the 

- 

minimum Q- seminorm is conditional on P least squares, so again X MX is 
X 

smaller for BLE. Notice that this estimate is very sensitive to errors 

in the Q- matrix as can be seen in the RSMPE it appears unsuitable 
X 

for our purposes. The other estimators are not affected by the 
non-positive definiteness of the model matrix. 

It should be noted here that a BLIMBE could be found for a posi- 
tive semidefinlte Q- matrix. In fact it has the same equation as for the 
positive definite case (3.3—41) (B. Schaffrin, 1982, private communica- 
tion) . It was tested for this 18-station experiment but proved to be 

quite unsatisfactory. The trace of the bias matrix was extremely small 
—8 2 

(~2 X 10 cm ) , so that one almost has an "unbiased" estimate. On the 
other hand, the standard errors of the parameter covariance matrix were 
unacceptably large with a magnitude of several meters and the deforma- 
tions were estimated very poorly. This is a classic example of not 
choosing an estimator by only its seemingly optimal properties without 
checking it also by simulations. Recall that the BLIMBE is formally 
minimum variance, but only in the class of minimum bias estimators. 

This is an extreme case of minimizing bias at the expense of variance. 
Therefore, a warning is issued for all potential users of the BLIMBE 

with a Q- seminorm. 

X 
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The good RSMPE results of the BLICUE are quite misleading as 
will become clearer in the zero deformation tests of section 4.3.5. 
Recall that 3 cm and 6 cm Gaussian noise were added to the "true" 
deformations in order to construct an incorrect model matrix The 

BLICUE constrains CX = CX (section 3.4.2) so that the random noise 
has not significantly changed the overall six CX values. Therefore, this 
is not the best way to test the BLICUE estimate. Later results will 
show that the constraints CX = CX will make the BLICUE estimate too 
sensitive to certain types of errors in the model matrix and therefore 
also unsuitable for our purposes. 

We are left then with the strong and weak Bayesian approaches, 
i.e., Bayesian versus BLE. It seems clear that the weak approach, that 
is introducing the expected deformations through the model moment 
matrix is preferable to the strong approach of direct use of the 
a priori information. When the model is "correct", the RSMPE 's of both 
estimates are at the same level (approximately 0.7 cm). The AMl-2 com- 
puted magnitudes of deformations for the two year period are listed in 
Table 8. They have an RMS of 6.1 cm so that most of the 
deformation is being recovered. On the other hand, the BLE is less 
sensitive to errors in the geophysical model. Note that for both esti- 

mates, the variance of unit weight is close to unity for the "cor- 

''2 

rect" model case. In the presence of model noise, cr^ indicates an 

imcompatibility between the geodetic observations and the geophysical 

model. One can then make a case that the Bayesian estimate may be more 

appropriate for testing a particular geophysical model, as indicated by 

''2 

the larger values for 0^. 
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4.3.4 MERIT-COTES Experiments 

The planned MERIT 83/84 main campaign (Wilkins, 1981) may be the 
first opportunity to establish the frame for a future CTS, considering 
that approximately 20 globally distributed stations will be available 
with a combination of the best VLBI, SLR and LLR instrumentation (Muel- 
ler et al., 1982). The strength of the reference frame is actually an 
indication of how well the polyhedron samples the earth. For a finite 
number of stations, then we investigate how "optimal" are the possible 
MERIT-COTES networks from the point of view of monitoring deformations. 

In (Mueller et al., 1982) it was assumed that no geophysical 
model is adopted and a free adjustment used for estimating deformations. 
The proposed MERIT-COTES (Working Group on the Establishment and Mainten- 
ance of a Conventional Terrestrial Reference ^stem) networks (see 
Tables 9, 10 and Fig. 5) were compared to optimal network (polyhedra) 
designs for different numbers of vertices. These optimal polyhedra re- 
sult from distributing p points on a sphere so that they are, in some 
sense, as far apart as possible from one another. It was stated there 
that an analogous optimality criterion is that the origin of the coordi- 
nate system defined by the p points is best determined (at least from 

the point of view of trace and determinant optimality defined below) . 

This is proven in Appendices B-D using some recent results in optimal 

design theory. It was shown that to a good approximation the above cri- 

terion also provides the best configuration for analysis of polyhedron 
deformations. Obviously, the distribution of stations is constrained 
by various factors, foremost of which is the location of the land masses. 


146 



147 


Fig. 5 MERIT-COTES Stations and AMl-2 Computed Velocities 
(see Table 10) 








Table 9 MERIT-COTES Global Networks 
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Table 10 HEBIT-COTES Stations and AHl-2 Velocities 
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The purpose of those simulations was to examine how close the possible 
networks could come to the Ideal case of being able to locate the sta- 
tions anywhere on earth. In other words, how close can one come to 
constructing an optimal polyhedra on the available land areas. 

In the case when a geophysical model for the station deforma- 
tions Is adopted, this must be a factor In the design of an optimal 
network. Unlike In the previous case, the optimal polyhedra are not so 
apparent. Here, we use as a basis for comparison the 18 station net- 
work distributed over the six major tectonic plates as described In 
ectlon 4.3.1. In addition, one station Is added to both the smaller 
Nazca and Arabian plates to form a 20 station network. In constructing 
these "optimal networks" we were guided by the considerations outlined 
In Section 4.3.1. The MERIT-COTES networks of (Mueller et al., 1982), 
now reanalyzed using the AMl-2 model, are compared to these two networks. 
We can then compare the different design measures under the two assump- 
tions, no plate model or AMl-2 model, to test if the conclusions in 
(Mueller et al., 1982) differ in any way. 

In (Mueller et al., 1982) the corresponding covariance matrices 
for the deformation estimate (3.3-55) and (3.3-83) are given assuming 
M = Q"^= I by 

= Oq n"*" (4.3-9) 

= Oq [(N + C^C)”^ - C^(CC^CC^)"^C] (4.3-10) 
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respectively, where 0^ = 1. Here we use the BLE for the comparisons 
based on the positive assessment of that estimator resulting from the 
simulations of the previous two sections, specifically the mean square 
error matrix 

MSE(X) = Qx - (3.3-147) 

A canonical notion of design optimality is not available, and 
we will briefly describe several common design criteria. All of these 
can be expressed in terms of the reduced elgenspace (the non-zero 
eigenvalues and their corresponding principal eigenvectors) of MSE(X). 

A-optimality is defined as minimizing the average variance 
(the A-measure) , in this case the average mean square error , or 
equivalently the spectral mean, i.e.. 


1 ^ 1 2 
min trace [MSE(X)] = min — 1 0 . 

3p dp 3. 


(4.3-11) 


3p - 6 

^i 


(4.3-12) 


where are the non-zero eigenvalues and 0^ the diagonal elements of 
MSE(X) . 

D-optimality is defined as minimizing the determinant of MSE(X) 
raised to the 1/ (3p - 6) power (the D-measure) or equivalently 

3p-6 ■ ^ 

n X 

i=l 


min 


3p - 6 


(4.3-13) 
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E-optimality is defined as 


min X 

X^ (4.3-14) 


We shall refer to the maximum eigenvalue as the E-measure. 

Another criterion, though used usually for determining the con- 
dition of a matrix rather than for optimal design, can be termed 
C-optimality, and is defined by 


min 


max 


min 


(4.3-15) 


the ratio of maximum and minimum eigenvalues. This is useful, since it 
is unitless and Independent of scale factor, i.e., baseline precision. 

It should be noted that all of these criteria are rotationally 
and translationally invariant (Isotropic and homogeneous) (Grafarend, 
1974). That is, only the relative distribution of the stations affects 
the optimal design. 

Besides these four optimal measures, we compute the "average" 

2 

bias, (4.3-8) squared to be in the same units (cm ) as the first three 
measures . 

In Table 11 the optimal measures of the MERIT-COTES 18A and 
20A networks are compared to their corresponding "optimal" 18 and 20 
station nets. Recall that each baseline is observed twice, once by 
VLBI and once by laser, or equivalently the baseline lengths are 
observed with 3/*^ cm accuracy. Note that both the 18A and 20A nets 
cover only 5 tectonic plates. As expected, the "optimal networks" 
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Table 11 Comparison of MERIT-COTES Networks Design Measures 
with 18-Station 6-Plate and 20-Station 8-Plate 
Networks 



Bxperiaent 



Design Beasures 




Mo. 

Oescciption 

A-Heasure 

D-Heasure 

B-fleasure 

C-Heasure 

Bias 


A1* 

A2** 

01* 

D2** 

El* 

E2** 

Cl* 

C2** 




(CB 

2 ) 

(CB^ 

) 

(ca 

2 ) 

(unitless) 

(ca^ ) 


HBBXT-COTBS 

with 

Single Type 

of 

Instruaent**** 






VLBX : 







7 

Extended Polaris 

16. S 

3.0 

6.2 

1.5 

153.9 

5.7 

120 

479 

2. 1 

10 

Exp. 7 * OSH 

3.6 

2.3 

3.2 

1.4 

15.3 

5.4 

16 

143 

1.7 

11 

Exp. 10 * Sao Paulo 

2.9 

1.8 

2. 8 

1.4 

10. 1 

4.8 

12 

1 16 

1. 1 

12 

Exp. 1 1 « Santiago 

2.6 

1.5 

2.5 

1.3 

8.6 

4.6 

11 

108 

1.0 




Laser : 







13 

Status in 1983 

2.5 

1.3 

2.3 

1.3 

10. 1 

4.4 

15 

38 

U.9 

14 

Exp. 13 * Santiago 

2.2 

1.2 

2. 1 

1.3 

7.1 

4.4 

11 

63 

0.8 


HEfiXT-COTES 

with 

Coabined Lasers 

and VLBI'S**** 



12 

Operational (9)*** 
Tt8X»SLB*LlB(aU) 


1.2 

1. 9 

1.3 

12.3 

4.7 

28 

51 

1.0 

18A 

T.o 

0.9 

0.9 

1.1 

4. 1 

4.2 

16 

38 

0.6 

18B 

18C 

S£i:!£inui 111 

4.9 

2.3 

1.1 
1. 1 

3.2 

1.7 

’•3 

1.2 

39.7 

16.3 

4.8 

4.6 

75 

43 

35 

33 

0.8 

0.7 

180 

20A 

Bxp^llvS2SLB (all) 

1.9 

0.8 

1.0 

0.8 

1.6 

0.8 

1.2 
1. 1 

6.8 

3.3 

4.5 

4.4 

19 

15 

33 

111 

0.7 

0.6 

208 

Priaarj (3) 

3.7 

1.1 

2.7 

1.2 

23.2 

4.9 

48 

115 

0.7 

20C 

Continental t9) 

2.1 

1.0 

1.6 

1.2 

12.9 

4.7 

36 

110 

0.7 

200 

Global (9) 

1.8 

1.0 

1.5 

1.2 

7.7 

4.7 

22 

117 

0.7 


Coaparison Hetworks 









18 

6- PLATE (ALL) 

0.8 

0.6 

0.8 

1.0 

1.5 

2.8 

9 

19 

0.4 

20 

8-PLATB (Abb) 

0.8 

0.6 

0.7 

1.0 

1.4 

2.9 

9 

19 

0.4 


* Ho plate eodel, H=X (3 ce baselioe accuracy) 

** AH 1-2 Hodel (3 ca baseliae accuracy — 2 year interval) 

•«* Huabecs in parentheses indicate nuaber of collocated sites 

***• for station locations see Table 9 and Figure 5 
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yield better results as reflected particularly in the E and C 
measures implying that the MERIT-COTES networks could be substantially 
improved. The large C-measure for 20A indicates that the two European 
stations added to 18A degrade the structure of the MSE matrix since 
these stations are added in a small region of a plate that is already 
densely covered in that particular area. Note that the bias term for 
the "optimal nets" is also smaller. 

It should be mentioned that the comparison of the no model and 
AMl-2 model experiments also shown side by side in Table 11 needs some 
qualification. First, two different estimators have been used. 

Second, the no model case is independent of the time interval between 
observations, while the AMl-2 Q matrix is a function of time. As can 

A 

be seen by examination of the results, the conclusions in (Mueller, 
et al, 1982) change in degree. That is, there are much smaller dif- 
ferences between the various collocation schemes, indicating that the 
absolute motion model provides an underlying frame of reference for 
the monitoring of deformations. 

Note that only the C-measures increase in the AMl-2 case. This 
is due to the positive semidefiniteness of Q as compared to the posi- 
tive definite M = I case. 

The conclusions for the single typed instrument experiments 
also do not change. Note though that the geophysical model improves 
the weaker VLSI experiment 7. This points to the advantages of a 
geophysical model. It is possible to perform a 3-dlmensional deforma- 
tion analysis in a network of less than global extent. For example, in 
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the VLBI experiment 7 stations are distributed over the North American 
plate and only the westernmost portion of the Eurasian plate. On the 
other hand, the free adjustment (M = I) requires a more global distri- 
bution for a meaningful 3-dimensional analysis. Furthermore, consider 
the polyhedron as wrapped completely about the earth. If the network is 
not closed, that is if there is a gap due to certain baselines not 
being measured, the free adjustment will have added to it further rank 
deficiencies of a geometrical nature besides the six that result from 
the coordinate system definition problem. This kind of singular situa- 
tion could occur for example in a limited VLBI network where problems of 
mutual radio source visibility make certain baselines nonobservable. 

A realistic geophysical model could be very helpful in overcoming such 
problems particularly using the BLE. 

The 12 station VLBI and the 13 (1^) station laser networks are 
again of basically the same quality. The structure of the laser net is 
stronger as reflected in the C-measure. This is most likely due to the 
distribution of stations over 5 plates Instead of 4 in the VLBI case as 
well as the better coverage of the plates on which the stations are 
located. 

Finally, note that the MERIT-COTES stations may not be in their 
optimal locations from the point of view of geophysically stable (on 
the Intraplate level) sites, for example, the Southern European and 
Japanese stations. Furthermore, it would be prudent, as mentioned 
before, to Increase the number of stations per major plate to 
four or five to increase the reliability of these subnets for interplate 
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motion detection. Intraplate and local motions could be monitored by 
filling in the subnets (within each plate) by, for example, GPS inter- 
ferometric observations (Counselman and Shapiro, 1978). 

4.3.5 Zero Deformation Tests 

In this section we test the situation in which the deformations 
are assumed to behave according to some model but actually there is no 
motion. This means that the difference between the baseline lengths at 
two epochs are due solely to observational noise. In these tests we 
use the Q matrix computed from AMl-2 but only 3 cm Gaussian noise is 
added to the initial baseline lengths to simulate the baseline observa- 
tions. We use the 18-station network. 

The results are listed in Table 12. As is apparent from the 
RSMPE, the biased estimates are less affected by the introduction of 
the faulty model, particularly the BLE. This is explained, as was 
noted before, from the deformations being introduced in a weak way in 
the biased estimates (MINDLESS and BLE) adjustments through the moment 
matrix. This same information is applied in a strong way as correc- 
tions to the station coordinates in the conditionally unbiased adjust- 
ments (Bayesian and BLICUE) . 

The results from this test point again to the BLE as the best 
deformation estimator. Recall that for the BLE estimation model, the 
deformation parameter vector X was assumed random and the model 
derived expected deformation X deterministic, while for the other 
three estimation models the opposite assumption was made. The BLE 
approach then is to filter the deformations (signals) from the 
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Table 12 Zero Deformation Simulations (&H1-2, 

3 cm baseline accuracy, 2 year time 
time interval, 18-station B-plate network) 


Estimate 

Property hiholbsS 

BLE 

Bayesian 

BLICUE 

yTpy 

(unitless) 

105.5 

147.5 

181.6 

104.6 

105-5 

X ^ MX 

(uuitless) 

11-7 

4.0 

81.0 

142.3 


93.4 

(1.3)e* 

28.7 

(0.7) 

0 

0 

j THJCOV(X))* 
<cm2) 

130.2 

(1.6)** 

- 

44.0 
(0-9) *♦ 

124.3 

(1-5)** 

TH jHS| ((X) ) ♦ 

223.6 

(1.9)** 

48.3 

(0.9)** 

44.0 
(0.9) ♦♦ 

124.3 

(1-5)** 

A2 

(unitless) 

t2 

(unrtless) 

1.0 

0-7 

1.0 

1.9 

1.0 

C-Measure 

(unitless) 

70.9 

19.4 

13.9 

20.0 

SHPE Cyar) 
(cm2) 

2.6 

0.5 

11.7 

23.2 

SilPE^ (^ias) 

0.0 

0.0 

0.0 

0.8 

BSHPE 

(cm) 

1.6 

0.7 

3.4 

4.9 


* A Priori values 

Numbers in parentheses are root mean square values (cm) 
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observational noise. This explains its success in this test since only 

noise was present in the simulated observations. The fact that the 

model matrix indicated expected deformations (with an EMS of about 

6 cm) did not alter the filtered deformations significantly. On the 

/v2 

other hand, the a posteriori variance of unit weight (Oq = 1) did not 
indicate anything peculiar although that there were no statistically 
significant deformations could have been determined from examining the 
deformation estimates. The Bayesian a posteriori variance of unit 
weight (Cq = 1.9) indicates that baseline observations and the geo- 
physical model were incompatible. This is seen again, then, to be the 
main positive point of the Bayesian estimator. 

4.3.6 No Model Tests 

In this section we consider the case where the AMl-2 model is 

correct but no model weight matrix is used, i.e., M = I. Alternatively, 

this could mean that we expect no secular deformations with a particular 

2 

uncertainty, or M = k I. Here we assume this to be at the 10 cm level. 

In this case, as described in Chapter 3, the BLICUE reduces to the 

BLIMBE, the Bayesian to the BLE which approximates the BLIMBE to a 

2 

degree depending on the scale factor k . Only the BLIMBE should be 
used in this case since the normal matrices can be quite ill condi- 
tioned for the other estimators as indicated by the C-measure in 
Table 13, yielding unstable covariance (or mean square error) matrices. 

The important result from this test is that under the assump- 
tions made (3 cm baseline noise, AMl-2 motions and a 2 year 
re-observation period) the RSMPE is at about the 5 cm level. 
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Table 13 No Plate flodel Siaulations (B=l) (3 ca baseline 
accuracTf 2 year tiae interval* 18-station 
6- plate network) 


Esti aate 


Property 

BLIHBE 

BLE 

Bayesian 

BLICUE 

Y^p? 

(unitless) 

118.8 

119.0 

119.0 

8.9 

118.8 

(unitless) 

9.2 

8.9 

8.9 

9.2 


600.0 
(3. 3) * ** ♦ 

601.7 
(3.3) ♦♦ 

0 

0 


84.5 
(1.3) ♦♦ 

- 

682.7 
(3.6) ♦♦ 

684.5 
(3. 6)*v 

TR j«S|^(X) ) ♦ 

684.5 
(3, 6) ♦♦ 

682.7 
(3.6) ♦♦ 

682.7 
(3.6) *♦ 

684.5 
(3. 6) ♦♦ 

A2 

(unitless) 

A2 

0 

(unitless) 

1.0 

0.2 

1. 1 

1.1 

1.1 

C- Measure 
(unitless) 

8.7 

201.0 

201.0 

200.0 

StfPE (Var) 
(cm2 ) 

22.7 

22.6 

22.6 

22.7 

SflPE^ (|ias) 

0.8 

0.8 

0.8 

0.8 

BSHPE 

4.8 

4.8 

4.8 

4.8 


(ca) 


* A Priori values 

** Nuabers in parentheses are root aean square values (ca) 


159 




Comparing this to the previous tests in which somewhat erroneous models 
were used, it can be concluded that it is preferable to apply a somewhat 
incorrect model than no model at all. However, this is contingent on 
the plate tectonic theory being essentially correct, though the precise 
motions may not be well known. 

4.4 The Treatment of Errors in 

In the simulations of this chapter, we assumed a linear, model 
for the computation of expected deformations from the station veloci- 
ties given by the plate tectonic model. Furthermore, we neglected the 
effect of uncertainties in the fundamental CTS coordinates X^. That 
is, we treated as errorless even though it will have been estimated 
from the observations of different measurement systems as described in 
section 2.2. Whether or not to consider the associated covariance 
matrix of X^,Z when computing deformations is a matter of philosophy. 

U Xq 

In any case, we seek a well defined datum to which deformations of the 
polyhedron are referred. It is clear that since the deformation prob- 
lem is dynamic, the Xq estimate can never be improved except by a rede- 
finition of the reference frame initial epoch at a later time when 
improved geodetic observational accuracy would warrant it. However, 
since we are primarily interested in the changes in the fundamental 
coordinates (X-X^), the effect of errors in X^ on these quantities 
should diminish with time. 

In this section, we add an offset term to the linear model 
(4.2-16) and indicate how errors in X^ can be incorporated into the 
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estimation by means of the model matrix This is done using a state 

transition matrix approach which is almost trivial for this case but can 
be generalized to deal with more complicated models. The motions of the 
polyhedron are referred to an initial epoch, i.e., to the fundamental 
polyhedron. If there were no deformations in the initial polyhedron, it 
would continue to rotate with the earth in its initial configuration. 
However, the earth is deforming so that perturbations are present which 
must be monitored. These perturbations may be periodic or secular. 

Consider that this motion could be described by a set of simul- 
taneous first-order differential non-linear differential equations 
(Llebelt, 1967) 

g = H(X(t),d(t),t) (4.4-1) 

where d(t) is a set of specified forcing functions. The integration of 
these equations of motion using the initial condition (state) X(tQ)=XQ, 
results in the deformed state (the state vector) X(t) at some later 
epoch. Since we invariably have to linearize our problems, we can 
simplify (4.4-1) by assuming a linear system 

g = F(t)X + G(t)u(t) (4.4-2) 

where u is a set of specified forcing functions which are related to 
the time rate of change of X by the matrix G which in our case would 
be computed from adopted earth models. In our application, we assume a 
homogeneous linear system so that u(t) = 0 and 


161 


# = F(t)x 


(4.4-3) 


This has a solution of the form 


X(t) = S(t,tQ)X(tQ> 


(4.4-4) 


where 


S(t,tQ) = F(t)S(t,tQ) 


(4.4-5) 


and 




(4.4-6) 


S is called the state transition matrix familiar in satellite orbit 
determination. If we can determine this matrix, then the deformed 
state of the polyhedron can be computed at any epoch by operating S 
on the initial state X^. In our case, we assume (compare to (4.2-16) 

X(t) = Xq + Xg(t-tQ) = Xq + V(t-tp) (4.4-7) 


where X^ denotes the fundamental coordinates and Xq the change in these 
coordinates computed from the AMl-2 plate tectonic model. The equa- 
tions of motion are easily written as 



(4.4-8) 
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Their solution is given in terms of S 


» 


and Xq by 
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0 
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• 

z 


0 

0 

0 

0 

0 

1 


• 

/o 


X 

= 



S(t 

1 *^ 0 ) 

• 



Xo(t) 


(4.4-9) 


Of course, this is just equivalent to (4.4-7). To compute the moment 
matrix Q 

A 

^ + (t-t_)^z* (4.4-10) 

and 

= E{XX'^} 

= E{(Xg + Xgt)(XQ + XQt)'^} (4.4-11) 

- E{XjXj} + (t-tj,)2E{i„xJ} 


where we assume that 


E{XqXJ} = E{XqXJ} = 0 (4.4-12) 

It follows that (compare to 4.2-17,18) 

T 

where the X^Xq term is already taken care of in the linearization of the 
deformation mathematical model (3.2-1). 
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In order to add the effects of errors In X^, then, to the 

deformation estimation. It Is necessary to add Its covariance matrix 

Z to the model matrix Q„. Thus, the linearization point X. of (3.2-2) 
Xq X U 

Is seen to be stochastic and therefore so Is L^, the Initial baseline 
lengths of the polyhedron and the zero order term of the Taylor's expan- 
sion. It follows that for (3.2-6) and (3.2-7) 

D[V] = = aQ(P“^ + P~^) = D[L] (4.4-14) 

4 2-1 

^ where <^qPq Is the covariance matrix of the Initial distances Lq pro- 
pagated from Z . The same dispersion matrix for V Is used for all 
the four estimators of Chapter 3. 

Although this model Is simple. It does provide a general method 
whereby a more complex model matrix Q„ could be derived. It could 
Include complicating factors such as Interplate and local motions, 
tidal effects, etc. This would require differential equations of the 
form (4.4-2). A finite element method may be most applicable to solving 
this system and determining the state transition matrix. In this case, 
the original state (X^) of the polyhedron would serve as Initial condi- 
tions and thete could also be appropriate boundary conditions. In this 
way, too, the deformations of non-CTS stations could be determined by 
densifying the finite element mesh. 
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5 CONCLUSIONS AND RECOMMENDATIONS 


5.1 Alternatives for Reference System Maintenance 

The frame of the future CTS Is to be defined at an Initial epoch 
by an adopted set of spatial coordinates of a global network of observa- 
tories mainly VLBI, SLR and LLR stations. In order to Insure that the 
earth orientation parameters are referred to the same set of axes, the 
deformations of the polyhedron need to be estimated periodically, l.e., 
a new set of CTS coordinates. This Is what Is referred to as main- 
taining the reference system so that the frame Is accessible to the 
user by the earth orientation (and translation) parameters In a con- 
sistent and accurate manner. 

In order to maintain the reference system on the deformable 
earth, either some constraints must be applied or geophysical models 
adopted, or both. One body of opinion holds that no geophysical model 
should be adopted at least not in the initial stages of the new CTS 
operations. During its early stages, one of the functions of the CTS 
could be to test If the geodetic observations are indicating motions 
compatible with those predicted by plate tectonic theory. Drewes 
(1982) presents an estimation procedure to estimate plate rotation 
parameters from a combination of geodetic and geophysical data, l.e., 
the inverse problem of Instantaneous plate kinematics (Minster and 
Jordan, 1974). Once the plate parameters are estimated, the forward 
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problem yields expected deformations of the stations that could be used 
to improve CTS deformation analysis. 

If no geophysical model is adopted, all the four estimators of 
Chapter 3 reduce to the I-norm BLIMBE, or to the familiar free adjust- 
ment. Therefore, one is limited essentially to applying the set of 
inner constraints CX = 0 which* without any real physical justification, 
impose no net rotation nor translation for the estimated deformations 
X . As indicated by the simulations of Chapter 4, if there are secular 
deformations of the CTS stations at the level predicted by the abso- 

/V 

lute motion plate models, the constraints CX = 0 are inappropriate. 

The reference system will nevertheless be maintained in a well defined 
manner. However, the CTS will then have fairly high sensitivity to 
changes in the distribution of the observing stations, and moreover to 
the actual locations of the chosen sites contradicting one of the 
requirements of Chapter 1. With time, distortions may accumulate in 
the system and the deformation estimates of the free adjustment may 
less and less resemble the physical deformations. 

Another body of opinion maintains that a geophysical model 
should be adopted from the initial stages of CTS operations. After all, 
one of the primary reasons for the establishment of a new CTS is the 
general acceptance of plate tectonic theory. This approach does not 
necessarily contradict the requirement of avoiding as much as possible 
dependence on geophysical hypotheses. The CTS frame is still defined 
by the coordinates of the fundamental polyhedron and is invariant with 
respect to an adopted geophysical model. The geophysical model does 
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affect the periodically updated CTS coordinates though hopefully 
Improving the estimation of station deformations as compared to the 
free adjustment. In order to reduce the dependence on geophysical hypo- 
theses, any deformation estimator should not be highly sensitive to the 
adopted model. Furthermore, if the model is incorrect, a good estima- 
tor will alert the investigator to this. 

Three alternatives have been suggested for maintaining the 
reference frame when a geophysical model is adopted and four estimators 

/s 

proposed. The first alternative uses the constraints CMX = 0 and its 
corresponding estimator is the BLIMBE, i.e., the weighted free adjust- 
ment. As shown in Chapter 2, the resulting reference frame axes are a 
discrete Tisserand's mean axes of crust. In the case of a positive semi- 
definite Qjj matrix (associated with any absolute motion plate model) the 
minimum bias estimator is unsatisfactory as mentioned in Chapter 4. 

The MINDLESS proposed in section 3. 3. 1.4 is found to be quite sensitive 

to errors in the Q matrix and, therefore, also unsatisfactory for defer- 
A. 

mation analysis. 

The second alternative is to combine the baseline measurements 
and the expected deformations, computed from an absolute motion plate 
model* without any constraints. This implicitely fixes the CTS frame 
axes in the mantle (and to the mantle and crust which rotate together in 
a mean sense). Two estimators presented follow this approach, the BLE and 
Bayesian estimators. The BLE introduces the geophysical model through a 
moment matrix in a weak Bayesian manner , while the Bayesian estimator does 
so directly in a strong Bayesian manner. When the model is correct, both 
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approaches provide good estimates. On the other hand, the BLE is less 
sensitive to the errors in the geophysical model. In the case of 
applying the AMl-2 model when there are in reality zero deformations 
(only observational noise) the BLE does not indicate any significant 
deformations. This follows from its being a filter rather than a true 
estimator. Note that this approach is less sensitive to the distribu- 
tion of the observing stations and in the frequency of observations 
since the geophysical model is good for any station location (although 
unstable areas should be avoided) . It would seem that conversely the 
dependence on the geophysical hypothesis should Increase. This is so 
for the Bayesian estimate but not for the BLE as described above. 
Furthermore, an extension of BLE (least squares collocation) can be 
used to predict deformations in case of trouble at a number of CTS 
stations, as well as to estimate the coordinates of new stations. 

/N _ 

The third alternative is to use constraints of the form CX = CX 

/s — 

(instead of CX=0) where X is computed from the geophysical model as 
explained in Chapter 4. It has been seen from examination of the cor- 
responding BLICUE estimator that this approach is a hybrid of the two 
previous ones. It has been shown to be quite sensitive to errors in 
the geophysical model in certain cases and does not give an indication 
of the presence of a poor model. This can be explained by its weighted 
constraint Interpretation. 

Summarizing, on the basis of the simulations of Chapter 4 it 
seems preferable to adopt even a weak but realistic absolute plate 
motion model for CTS operations than none at all. In this case, the 
BLE seems most suitable for deformation analysis. If the model is of 
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poor quality, the estimation algorithm should indicate this. If no 
model is adopted, the I-norm BLIMBE can insure a well defined reference 
system but of questionable physical significance. ^ 

5.2 MERIT-COTES Networks and the 18-Station 6-Plate Network 


Several conclusions can be drawn from the comparison of the 
MERIT-COTES networks to the 18— station six -plate network of section 4.3.4 
other than the ones already mentioned. The planned MERIT-COTES net- 
works for the MERIT main campaign provides a good starting point for 
the eventual establishment of a new CTS. However, the distribution of 
stations over the major plates should be improved. Furthermore, the 
number of stations per plate should be at least four from the point of 
view of reliability as well as geometrical strength (for which three 
seems to be adequate though). Nevertheless, it is strongly recommended 
that over the duration of the MERIT main campaign two or three short 
sessions be devoted to observations from all stations. The first such 
session at the start of the campaign could be used to establish an 
initial CTS frame. Subsequent sessions could monitor deformations 
using the recommended BLE estimator. For testing the compatibility of 
the baseline measurements with geophysical models, particularly the 
plate motion hypotheses, the Bayesian estimator is more appropriate. 

Under the following assumptions, 

1. an 18-station six-plate network 

2. 3 cm baseline length accuracy (this implies removal of the 
larger systematic errors) 
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3. a fairly reliable absolute motion model 

4. no significant anomalous deformations 

It can be concluded from the simulations of Chapter 4 that station 
deformations could be estimated at the 1 cm level from periodic reobser- 
vation of baseline lengths every one or two years. This will essen- 
tially remove the effects of Interplate motions from the error budget 
of the required 5 cm accuracy short term variations In polar motion and 
earth rotation. The resulting reference frame will provide an accurate, 
well defined and consistent zero-order global network for geodetic and 
geophysical studies. 

A potentially major problem In this optimistic scenario Is In 
assumption 4 above. This could be due to site stability problems. 
Intraplate motions and unmodeled tidal effects such as ocean loading. 
Therefore, the sites must be chosen very carefully by geophysical sur- 
veys. Once they are chosen, local effects could be monitored by on-slte 
observations such as gravity observations and local geodetic surveys. 
Intraplate motions and some local effects could be monitored by GPS 
derived baselines, particularly In the Interferometric mode. Further 
Investigations are recommended to study optimal ways to Incorporate 
these observations In a well defined manner Into CTS operations. 

5.3 A Final Comment on Estimation 

The examination of the four estimators In Chapter 3 and the 
simulations of Chapter 4 have provided Insight Into the meaning of 
biased and unbiased estimation. An estimate of a parameter vector Is 
unbiased If and only If there are no restrictions on the estimation. 
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In order words, the estimate can conceivably take on any value In the 
parameter space. 

If the design matrix A Is singular, there exist relationships 
among the parameters, restricting the resulting estimate. In the pres- 
ence of a priori information on the parameters, it is still possible to 
construct an unbiased estimate (strong Bayesian) with the assumption 
that the expectation of the a priori parameter vector is equal to the 
true one. (Under this assumption, it is unbiased also when A is of 
full rank.) Formally, we have constructed an unbiased estimate but we 
have "biased" the estimate in the direction of the a priori value. 
Strictly, this estimate can be only considered conditionally unbiased. 

If the assumption is correct (within its uncertainties) , the result will 
be an estimate with lower variance than the BLUE (when A is full 
rank) or BLIMBE (when A is rank deficient). If incorrect, the esti- 
mate will suffer accordingly. This is the basic danger in applying the 
strong Bayesian estimate. 

A better approach, in the case of a priori information and 

whether or not A is of full rank, is to use a biased estimate (BLE) , 

by constructing an appropriate moment matrix Q . This reduces the 

A 

effects of an incorrect assumption for X , which follows from its fil- 
tering interpretation (treating X as a random variable) . Its more 
general prediction capabilities are also very useful. 

Alternatively, in the case of singular A , the BLIMBE provides 
a minimum bias estimate (in contrast to minimum mean square error for 
BLE) that might be preferable depending on the application. For 
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non-singular A , the BLIMBE reduces to the BLUE invariant with respect 
to a chosen weighted norm in the parameter space. 

Summarizing, if even weak but essentially realistic a priori 
information on the parameters is available, it makes little sense to 
ignore it. For non-singular A , this is what we do though when the 
the BLUE is used. However, if we do not choose to ignore this informa- 
tion, why try to construct a seemingly unbiased estimate (the Bayesian 
estimate) . Perhaps the term unbiasedness (and the related estimability) 
is conveying Incorrect connotations of approbation that are difficult 
to shed. In the case when the BLUE does not exist (as in the deforma- 
tion problem) we may be tempted to apply a conditionally unbiased esti- 
mate (Bayesian again or BLICUE) in order to seemingly make up for the 
biasedness (or non-estlmability) of the problem. We have shown that 
(at least for analyzing deformations) it is better to stay within the 
class of biased estimators in the application of a priori information. 
Perhaps this conclusion holds for problems that are non-singular to 
begin with. That is, in the presence of a priori information, it would 
be a better practice to move into the class of biased estimators. This 
seems to conform with the current mode of thought and investigations of 
the broader statistical community. (See (Trenkler, 1981) for a good 
review. ) 
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APPENDIX A 


THE WEIGHTED PSEUDOINVERSE 


A. 1 Conditions for a Weighted Pseudoinverse 

Consider an Inconsistent set of linear equations Y = AX, for a 
rank deficient A matrix. A Is of dimension n x u and maps X In E^ 
Into Y In e’^. A weighted (ellipsoidal) norm Is defined In E^ 

II X ||„ . (xV^''^ 

and In e’^ 

II L lip = (L^PL)^/^ 

We assume that P and M are positive definite matrices. The fol- 
lowing theorem from (Rao and Mltra, 1971) provides the conditions for a 
solution X = GY to be minimum M-norm P least squares. 

Theorem . Let there exist a matrix G such that GY Is a minimum M-norm 
P least squares solution of AX = Y. Then It Is necessary and sufficient 
that the following conditions hold 


AGA = A 


GAG = G 
(GA)"^M = MGA 
(AG)'^P = PAG 


(A-1) 

(A-2) 

(A-3) 

(A-4) 
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where P and M are positive definite matrices. 

A matrix G that fulfills these conditions is denoted by ^ mini- 

mum M-norm P least squares g- inverse or a weighted pseudoinverse (Bouillon 
and Odell, 1971). In the case M and P=I this reduces to a"*^, a minimum 
norm least squares g-inverse, or simply the familiar pseudoinverse. 

Consider conditions (A-1) - (A-4) when M = P = I. If a matrix 
fulfills (A-1), it is called a generalized Inverse A ; (A-1) and (A-2), 
a reflexive Inverse A ; (A-1) - (A-3) , a left weak inverse A ; (A-1) , 
(A-2) and (A-4), a right weak inverse A ; and all four conditions, a 
pseudoinverse. Only a pseudoinverse is unique. These relationships are 
illustrated in Fig. 6. ' 

Conditions (A-1) and (A-4) are equivalent to 

A^PAG = a’’p (A-5) 


which can be proven as follows. Assume that (A-5) is fulfilled. Then, 
(AG)'^P = g'^a'^P 

rp m T T* T* 

= G A PAG = [G A PAG] 

= [G^a'^P]'^ = PAG 

which is just condition (A-4). Using this result 

AGA = P"^PAGA 

= P”^(AG)^PA = p'^G^a'^PA 

= [A^PAGP"^]"^ 

= [a’^pp"^]'^ 


= A 


which, is condition (A-1) , 
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Conversely, given (A-1) and (A-4) 


T T T 

A PAG = A^(AG)^P 

= A G A P 


(PAGA) 
T 


= (PA) 
T 

= A^P 


which completes the proof. 

Condition (A-2) and (A-3) are equivalent to 


T T 

G MGA = G M 


The proof is as follows. Given (A-6) 


T T T 

(GA) M *= A G M 

T T 

= AG MGA 

= (a^g'^mga)'^ 

= (aVm)'" 


= MGA 


which is condition (A-3) . Using this result 

GAG = m"^GAG 

-1 T -ITT 

= M (GA) MG = M A G MG 

T -IT 
= [G MGAM ] 

T -1 T 
= [G MM ] 

= G 


(A-6) 
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which is condition (A-2) . Given (A-2) and (A-3) 

g\ga = g'^(ga)^k 

= [MGAG]"^ 

= [MG]*^ 


which completes the proof. 


A. 2 A Proof for A^^^ (M Positive Definite) 


Here we prove that G = M ^N(NM ^N)®A^P of section 3. 3. 1.1 and 
3. 3, 1.2 is A^^* For this we will need the results (Rao and Mitra, 
1971, p. 22) 

(a) One choice of (A*^)® is (A®)^ 

(b) A(A^PA)®A^PA = A and (a'^^PA) (a'^PA) ®A^ = A^ for any matrix P 

T 

such that R(A PA) = R(A) which automatically holds if P is 
positive definite. 

First, we show that G fulfills condition (A-5) , 
a'^PAG = NG 

= nm~^n(nm~^n)®a'^p 
= nm"^n(nm“^n)®nn®a’^p 

K T 

= NN®A P 
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For condition (A-6) 


G^MGA = PA[(NM"^)®]'^NM~^^”^N(NWN)®N 

= pa[nm”^n]%m“^n[nmn]% 

= PA[NM"^N]% 

= pa[nm“^n]®nm“^ 

T 

= G M 

Therefore, it follows that G is a minimum M-norm P least squares 
g-inverse. 

A. 3 An Interesting Relationship 

Here we show that G ® M ^N(N is an for the system of 

IM 

consistent normal equations NX = U. In order to prove this we show that 

conditions (A-1) - (A-4) are fulfilled (in this case P does not appear 

T 

in the conditions, being contained in N = A PA) using again (a) and (b) 
of A. 2. 

(1) NGN = NM"^N(NM"^N)^ = N 

(2) GNG = m"^N(NM“^N)®NM“^N(NM”^N)® 

= M~^N(NM"^N)® = G 
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(3) (GN)'^M = (MGN)'*^ 

= [N(NM“^N)%]'^ 

= N[NM~^)% 

= mm“^n(nm"^n)% 

= MGN 

(4) (NG)'^ = [NM~^N(NM“^)®]’^ 

= NM"^(NM~^N)® 

= NG 


proving the assertion for G . This proves that 


+ + T 


A. 4 A Proof for (M Positive Semidefinite) 


Here we prove that 


G = (N+M)“^N[N(N+M)“^N] ®A^P 


(A-7) 


of (3.3-119) fulfills the four conditions of section 3. 3. 1.4. It follows 
that G = A^j^ is, in this case, a minimum M-semlnorm P-least squares 
solution for AX = Y. We assume that P is positive definite and use 
(a) and (b) of A. 2. Then, 
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(1) PAGA = PAN%GA 

= PAN%(N+M)”^NIN(1W-M)“^N]% 

= PAN%A 
= PA 

(2) MGAG = M(N+M)"^N[N(N4M)”^ ]®N(N4M)"^N[N(N4W)"^N]®a'^P 

= M(N+M)“^N[N(1H^I)"^N]®a'^P 
= MG 

(3) MGA = (N+M)GA - NGA 

= N[N(N4M)”^N]%- N(N+M)”^N[N(N+M)"^N]% 

= N[N(N+M) ^]% - N (symmetric) 

= [N(N+M)“^1% - N]*^ 

= (MGA)’’ = (GA)’’^M 

(4) PAG = PA(N+M)“’'N[N(N4M)”’'n]®a’^P 

= PA(N4M)~’'N[N(N+M)“’’N]^NN®a'^P 

= pan%n®a’’p 

= PAN^a’^P (symmetric) 

= [pan^a’^p]’^ 

= (PAG)’^ = (AG)^^P 
which completes the proof. 
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APPENDIX B 


APPROXIMATE THEORY FOR OPTIMAL DESIGN 

2 -1 

Consider the Gauss-Markof f model (LjAXjOqP ) for the observation 
equations 


L = AX + V (B-1) 

where the rank of A is full. Each row of A is a u-dimensional 
vector in the design space, a subspace of e'^. The normal matrix N is 
known as the Fisher Information matrix (Federov, 1972) . Actually the- 
definitlon of Fisher information is much more general (see for example 
Silvey, 1980) but simplifies to N for the linear estimation problem 
given above. The design problem is choosing n-vectors A^,...,A^ 

Such that the covariance matrix (without loss of generality we assume 
P = I) 


n - 

= [ 2 A,A,]“-^ (B-2) 

^ i=l ^ ^ 

is minimized in some sense. Alternatively, we wish to minimize the 
information matrix N . Several design measures have been presented in 
section 4.3.4. For the linear model, the design matrix is independent 
of X . However, for the linearized model, it is dependent on 
the approximate values of X. In the case of deformation 
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analysis, this value Is given by Xq, the fundamental coordinates. 

Since the expected deformations are small, several cm/year over base- 
lines of Intercontinental extent, the theory developed for linear models 
is applicable here, too. 

Optimal design criteria, in general, consist of different func- 
tions of the information matrix 

(KN(d)) (B-3) 

where d denotes the vectors that make up a particular design. An 
optimal design, denoted by d , will maximize (j). Since d consists of 
a discrete number of vectors A^, it is not practical to apply standard 
optimization techniques to maximize (j). A way to generalize this prob- 
lem in terms of continuous functions has been developed by (Kiefer, 

1974) which he has termed approximate theory. We outline this approach 
as presented by (Silvey, 1980). 

Consider a design space , a compact subset of e'^. Only cer- 
tain collections of A^^eA can be considered as valid designs (for 
example, the number of vectors must exceed the number of parameters) . 

The collections or events can have probabilities assigned to them. The 
class of possible events forms a field, particularly a Borel field since 
any combination of these events is also an event and belongs to the 
field. Let H denote the class of probability distributions over the 
Borel field of possible events. That is, each possible design has a 
probability associated with it. Then ri£ Hcan be thought of as a 
design measure. For every n define its information matrix 
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V 


N(n) = E {A^A^} (B-4) 

where A^<£"A is a random vector with distribution n- The set of all 
possible information matrices is given by 

M = ^N(n) ; n e H} (b-5) 

Define (|) as a real-valued function over the set of u x u S 3 nnmetric 
non-negative definite matrices. The design problem can then be stated 
as determining Ti that maximizes <J>{N(n)} over H . Such a design is 
called 4)-optlmal. The properties of (|>, are outlined in (Federov, 1972; 
Silvey, 1980). The above definitions, then, allow us to consider con- 
tinuous optimal designs. 

Two directional derivatives are defined. The Gateaux deriva- 
tive of <() at N^, in the direction of N 2 is 

G^(N^,N2) = lim 7 { 4 >(Nj^+ N2) “ <|>(N^)} (B-6) 

e-K)'*’ 

If is differentiable at 

G^(N^,Za^N^) = Ia^G^(Nj^,Nj^) (B-7) 

a property that slmplfles the design problem considerably. The 
Frechet derivative of <}) at N^, in the direction of N 2 is 

lim ^ {4)[(1 -e)N3^ + eN2] " (B-8) 

By definition. 
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F^ttl,N2) - G^(Sj.H2-N2) 


(B-9) 


Two theorems are important in deciding whether a particular design is 
optimum. 

Theorem 1 . When (j) is concave onj^, n is (|)-optimal if and only if 

Fjj,{N(n*) , N(n)} < 0 for all nC H 

/ 

ic is 

Theorem 2 . If 4> is concave on /j and differentiable at n(H ) then r| 
is (j)-optimal if and only if 

F^{N(n*), A^A^} < 0 for all A^G A 

For similar theorems treating the more general case of a rank deficient 
A matrix, see (Silvey, 1980). Theorem 2 will be useful in the study of 
optimal polyhedra on a sphere (Appendix D) . In order to apply this 
theorem to this problem it will be necessary to compute the Frechet 
derivative for each <{)-optimality criterion. We do this for D- and 
A-optimality. 

For D-optlmality we define 

<J) = log[det N] (B-10) 

to ensure concavity. Then 

J = log det + eN^) - log det 
= log [det(Nj^+ eN 2 >N^^] 

= log [det(I+ eN 2 N~^)] 
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Expanding in a McLauren series 


J = log {detl+ ^ [det (1+ eN 2 N"^)]}^^q + 6(e^) 

but 

(det A) = det A tr (A ^ 

(Bodewig, 1959). Therefore, 

J = log{det 1+ e [det(I + eN 2 N~^) 

•tr((I+ ^ 

= log{l + e[det(I + £ N 2 N"^) tr((I + e N 2 N]^^)"^N 2 N];^) ] 

= log{l + e[detItr(N 2 N*^)]} + o(e^) 

= log{l+e tr(N 2 Nj^^)}+ o(e^) 

Since 

log(l + X) = X - J X^+... 

J =E tr(N2N]^^)+ o(e^) 

and we can compute (B-6) for D-optimality 

G (N^,N 2 ) = lim J etr(N 2 N^^) 
e^'*’ 

= tr(N2N^^ 

From (B-9) 


(B-11) 
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F^(N^,N2> = tr(N2-Nj^)N 




= tr[N„N7^- I ] 
2 1 u u 

= tr(N 2 N"^) - u 


(B-12) 


and 


'•/"i-Vi) 


tr(A^A^Nj^^) 

tr(A^N]^^A^) 
-1 T 

A±Ni A^ - u 


- u 


- u 


(B-13) 


where A^ is a row of the design matrix A . Then according to the 

It 

second theorem above, n Is D-optlmal If and only If 


A^N"^(n*)A^ < u for all A^G Si 


(B-14) 


For A-optlmallty we consider the concave function 


<()(N)= -tr(N”^) 


Then, 

K = <j)(N^+ EN 2 ) - (t>(Nj^) 

= -tr[(N^+ eN 2 )”^I + trN~^ 

= trN^^- [tr(N^+ eN 2)~^]^^0 + trN^^ + o (c 

Now, 

^ trA - tr ^ A 

and 
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= _A "1 M -1 
^ 3 e ^ 


(Federov, 1972). Then 

K = etr[(N^+ eN 2 )”^N 2 (N^+ £^ 2 )"^]^^^ + o(e^) 
= tr[N^^N2Nj^^] +o(e^) 


from which (B-6) for A-optimality is 


G^CNf.Na) = lim ^etr(N~^N 2 N]^^) 
e-K)"*" 

= tr(N~^N 2 N“^) 


(B-15) 


and from (B-9) 


F(Jj(N^,N 2)= tr[N^^(N2~Nj^)N]^^] 

= tr[N"^N2N^^-N^^] (B-16) 

- 

It follows, 


F^(N^,a][a) = tr(N]^^A^A^N"^) - tr(N”^) 
= tr(A^A^N~^”^) - tr(N^^) 

= tr(A^N"^"^A^) - tr(N"^) 

-1 -1 T -1 
= A^N^ Nj^ A^ - trNj^ 


(B-17) 


Therefore, a necessary and sufficient condition that a design ri is 
(j)-optimal is that 

A^N”^(n*)N“^(n*)A^ £ trN"\Ti*) for all A^G fi (B-18) 
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APPENDIX C 


POLYHEDRA 


C.l Some Polyhedra Definitions and Classifications 

A polyhedra is a finite set of polygons arranged in space in 
such a way that every side of each polygon belongs to just one further 
polygon, with the restriction that no subset has the same property. 

The polygons are called faces, the sides edges, and the juncture of 
several edges, vertices. In order to define the regularity of a poly- 
hedron, the concept of vertex figure is needed. A vertex figure of a poly- 
gon is the segment joining the midpoints of any two adjacent edges. The 
vertex figure of a polyhedron is the polygon whose edges are the vertex 
figures of all the faces that surround a vertex, i.e. the polygon formed 
by joining the midpoints of the edges which meet at a common vertex. 
Generally, this is a skew polygon. A regular polyhedron is one whose 
faces and vertex figures are all regular. Such polyhedra have faces of 
all of one kind of congruent regular polygons. (See (Coxeter, 1963; 

Fejes Toth, 1964) for more details.) 

The following scheme is used to classify regular polyhedra 

{p.q } (C-1) 

where p is the number of edges on each face and q, the number of vertex 
figures (also, the number of edges meeting at each vertex) (Coxeter , 1963). 
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The only possible (p,q)'s are for the five Platonic solids 

(3.3) - tetrahedron 

(3.4) - octahedron 

(4.3) - cube 

(3.5) - icosahedron 

(5.3) - dodecahedron 

See Fig. 7. All (p,q)’s can be inscribed in a sphere. They are also 
reciprocal (dual). This means that if we join the centers of adjacent 
edges by segments, we obtain the reciprocal polyhedron (q,p). This fact 
was used to construct near-optimal polyhedra (called dual polyhedra) in 
(Mueller, et al, 1982). 

A semi-regular polyhedron has regular polygons as faces , but the 
faces are not all of the same kind. There are 13 such polyhedra called 
after Archimedes. These all denoted symbolically by the number of the 
edges about one vertex. For example, for a polyhedron of 24 vertices 
we have (using the classification of Fejes Toth (1964)) 

(3. 3. 3. 3. 4) - snub cube (snub cuboctahedron) 

i.e., about each vertex we have 4 triangles and one square (Fig. 7) 
and 

(3,8,8) - truncated cube 

i.e., around each vertex there is 1 triangle and 2 octagons. For 8 
vertices, the antiprism (Fig. 7) is represented by 

(3. 3. 3. 4) 
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20 * 



Dodecahedron 


24 



Snub Cube 
Fig . 7 Examples 



Antiprism 


14 * 




Dodecahedron + Icosahedron 


♦ near optimal 

•♦minimum distance maximized - 
Fejes Toth , 1964 .Regular Figures 


of Optimal** Polyhcdra 
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i.e., about each vertex we have three triangles and one square. All 
three of these semlregular polyhedra can be inscribed in a sphere. 

C. 2 The Golden Proportion (Section) and Polyhedra Coordinates 


The pentagram is a five-pointed star constructed by extending the 
edges of a regular pentagon. The 10 outside edges of a pentagram is T 
times the length of the edges of the original pentagon, where 


+ 1 

— 


(C-2) 


is called the golden proportion. It is a root of 

- X - 1 = 0 (C-3) 

Many interesting relationships are attributed to t (Pugh, 1976) . Some 
of these relate to polyhedra. Dividing each edge of an octahedron by 
T : 1 yields an icosahedron. If a diagonal is drawn across each penta- 
gonal face of a dodecahedron, a cube is formed. The edges of the cube 
will be T times the edges of the dodecahedron. 

The Cartesian coordinates for the regular and semi-regular 
polyhedron described in the previous section are often expressed in 
terms of T. For the tetrahedron (p = 4) , a set of coordinates is 

(1,1,1) , (1,-1,1) , (-1,1,-1) , (-1,-1,!) (C-4) 

those of the octahedron (p = 6) 

(±1,0,0) , (0,±1,0) , (0,0, ±1) (C-5) 
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and those of the cube (p = 8) 


(± 1 ,± 1 ,± 1 ) 


(C-6) 


(±a,0,±h) , (0,±a,±h) 


/2 V, 1 

= ; h = — 


(C-7) 


For the icosahedron (p = 12) 


(0,±T,±1) , (±1,0,±T) , (±T,±1,0) 


(C-8) 


and the dodecahedron (p = 20) 


(0,±t'’^,±t) , (±T,0,±t“^) , (±t"^,±T,0) , (±1,±1,±1) (C-9) 


(Coxeter, 1963). For the semi-regular polyhedra, the antiprism (p = 8) 
Inscribed in the unit sphere has coordinates 


(±a,0,h) , (0,±a,h) , (±— ,±— ,-h) 


(C-10) 


where 


2 , 1 
— ■ - ; h = ■ - 

/ 1 + 2 / 2 ' 


(C-11) 


The snub cube has coordinates (w. McWorter, 1982, private communication) 


(±a,±l,±a^) , (±l,±a,±a^) , (±a^,±a,±l) 


(C-12) 


where =.5436890127... is a root of 


X^+X^+X-1=0 
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and those of the truncated cube 


(±a,±b,±b) , (±b,±a,±b) , (±b,±b,±a) 

where 



(C-14) 


(C-15) 
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APPENDIX D 


POLYHEDRA DESIGN 

"On a planet, say, ten Inimical dictators govern. How must the 
residences of these gentlemen be placed In order to be as far 
as possible from one another?" (Fejes Toth, 1964). 

Consider the following variation of the reference frame problem. 
A total of p polyhedron stations can be distributed over the earth's 
surface. Assume a spherical earth and that the stations can be located 
anywhere on the surface. We define the optimal reference frame as the 
one that best defines the center of the sphere-, l.e. the origin. Assume 
that the "observations" are the radii R to each polyhedron station. 

Our mathematical model Is then 

R^ = [X^-Xq)^ + (Yj^-Yq)^+ (D-1) 

The parameters are the center of the sphere (X^jY^jZ^). The optimal 
design Is then the choice of p stations, l.e. (X^,Y^,Z^, 1 = 1, p) . 

This Is analogous to the problem of how to distribute the stations that 
will best monitor the change In the size and shape (deformation) of the 
polyhedron. In both problems the orientation of the reference frame 
axes Is arbitrary. 

Linearizing the above model yields the elements of one row of 
A of' the design matrix 
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(D-2) 


. Pi-^0 "i-M 

“R— ~~T-j 

Assuming a unit sphere and that the origin is at (0,0,0) 


*1 ' '’‘i ifi "i’ 


a>-3) 


The information matrix N for this problem for a particular 
design, r), and assigning equal probabilities to each point, is 


N(n) = 


z X A a;: 

i=l ^ 


p 




Sym 



ZX^ Z_, 1 

i i 

i 1 

ZY 2 

ZY^X^ 

i 

i 1 

ZZ^2 


(D-4) 


As an example, consider the case of p = 8 and check the cube and anti- 
prism configurations for D-optimality and A-optimality. One set of 
coordinates for the cube or the unit sphere is given by the 8 combina- 
tions (C-6) 




We take this as a possible design t]q whereby 


= 24 


8 0 0 
0 8 0 
0 0 8 




Now, for Oq to be D-optimal (B-14) 


A^N"^(t1o)aJ < 3 
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for all . But is restricted to points on the unit sphere. So 

A^n"^A = [X^ Z^]3I 

2 2 2 

= 3(X^ + Y^+Zp £ 3 

indicating that the cube distribution is D-optimal. 

For A-optimality (B-18) 

A^N“^(riQ)N"^(nQ)A^ < tr (N"^(n)) for all A^ G /I 

We have 

= 9(X^ + Y^ + zJ) £ 9 

so that the cube distribution is also A-optimal. 

The antiprism can be shown to maximize the shortest distance 
between any two vertices of the polyhedron (Fejes Toth, 1964) . We test 
whether this criterion passes the D- or A-optimality test. One possible 
set of coordinates of the antiprism is given above by (C-10) . Then 


tX^ Zp91 


^i 


lZ.V 
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Is this less than 3 for any point on the unit sphere? Try the point 
(0,0,1). In this case 


^ = 3.82 > 3 
h 

which shows that the antiprism is not D-optimal. A similar calculation 
shows that it is not A-optimal, either. These two results are verified 
in the experiments of (Mueller, et al, 1982). However, the antiprism 
gave a lower E- and C-measure (Section 4.3.4) than the cube. 

Consider the dodecahedron distribution for p = 20. The informa- 
tion matrix is given by 
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0 0 0 


0 0 

N(no> - } 

1 

0 2 0 
T 

4 

0 0 0 


2 

0 0 T 


o 

o 

1 


/ 

1 

0 

“ 

0 


1 

0 

0 


CM 

+ 

CM 



0 

0 


0 

2 

T 

0 

4 

0 

1 

0 

4 

0 

T 

0 


0 

0 

0 


0 

0 

1 


0 

0 

T^+^+2 






““ 




Lm, 


- 


12 1 

j + ^+ 2)1 

T 


for R=v^ and therefore. 


NCtIq) = ^ (t^ + -^ + 2)1 


for the unit sphere, 


2,1 

^ ■'"z 

T 


But 

3 


so 

NCtIq) = yl 

and 

n"^(t1q) = 3 1 

which proves that this distribution is D-optimal. In fact, it has the 
same information matrix as for the cube. Therefore, A-optimality holds, 
too. It can be shown that all regular polyhedra (i.e. the Platonic 
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solids) have the same information matrix and therefore all are A- and 
D-optimal. This also holds for the dual combinations 

p = 14 : cube and octahedron 
p = 32 : icosahedron and dodecahedron 

(see Fig. 7) 
and for 


p = 18 : octahedron and icosahedron 

Thus, regularity (or semi- regularity) is the criteria for A- and 
D-optlmallty for the origin definition problem. 
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